
!------------------------------------------------------------------------!
!  The Community Multiscale Air Quality (CMAQ) system software is in     !
!  continuous development by various groups and is based on information  !
!  from these groups: Federal Government employees, contractors working  !
!  within a United States Government contract, and non-Federal sources   !
!  including research institutions.  These groups give the Government    !
!  permission to use, prepare derivative works of, and distribute copies !
!  of their work in the CMAQ system to the public and to permit others   !
!  to do so.  The United States Environmental Protection Agency          !
!  therefore grants similar permission to use the CMAQ system software,  !
!  but users are requested to provide copies of derivative works or      !
!  products designed to operate in the CMAQ system to the United States  !
!  Government without restrictions as to use by others.  Software        !
!  that is used with the CMAQ system but distributed under the GNU       !
!  General Public License or the GNU Lesser General Public License is    !
!  subject to their copyright restrictions.                              !
!------------------------------------------------------------------------!


C RCS file, release, date & time of last delta, author, state, [and locker]
C $Header: /project/yoj/arc/CCTM/src/aero/aero5/coags.f,v 1.4 2011/10/21 16:10:13 yoj Exp $

C:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
      SUBROUTINE GETCOAGS( LAMDA, KFMATAC, KFMAT, KFMAC, KNC, 
     &                     DGATK, DGACC, SGATK, SGACC, XXLSGAT, XXLSGAC, 
     &                     QS11, QN11, QS22, QN22, 
     &                     QS12, QS21, QN12, QV12 )

c  Calculates the coagulation rates using a new approximate algorithm
c  for the 2nd moment.  The 0th and 3rd moments are done by analytic
c  expressions from Whitby et al. (1991).  The correction factors are also
c  similar to those from Whitby et al. (1991), but are derived from the
c  Gauss-Hermite numerical quadratures used by Binkowski and Roselle (2003).
 
c  Revision History:
c   FSB 08/25/03 Coded by Dr. Francis S. Binkowksi
c   FSB 08/25/04 Added in-line documentation
 
c  References:
c   1. Whitby, E. R., P. H. McMurry, U. Shankar, and F. S. Binkowski,
c   Modal Aerosol Dynamics Modeling, Rep. 600/3-91/020, Atmospheric 
c   Research and Exposure Assessment Laboratory, 
c   U.S. Environmental Protection Agency, Research Triangle Park, N.C., 
c   (NTIS PB91-161729/AS), 1991
c
c   2. Binkowski, F.S. an U. Shankar, The Regional Particulate Matter
c   Model 1. Model decsription and preliminary Results, Journal of
c   Geophysical Research, 100, D12, pp 26,191-26,209, 
c   December 20, 1995.
c
c   3. Binkowski, F.S. and S.J. Roselle, Models-3 Community 
c      Multiscale Air Quality (CMAQ) model aerosol component 1:
c      Model Description.  J. Geophys. Res., Vol 108, No D6, 4183
c      doi:10.1029/2001JD001409, 2003.
C-----------------------------------------------------------------------

      IMPLICIT NONE

C *** Arguments:
      
      REAL( 8 ) :: LAMDA   ! mean free path [ m ]

                           ! coefficients for Free Molecular regime:
      REAL( 8 ) :: KFMAT      ! Aitken mode
      REAL( 8 ) :: KFMAC      ! accumulation mode 
      REAL( 8 ) :: KFMATAC    ! Aitken to accumulation mode

      REAL( 8 ) :: KNC     ! coefficient for Near Continuum regime

                           ! modal geometric mean diameters: [ m ]
      REAL( 8 ) :: DGATK      ! Aitken mode
      REAL( 8 ) :: DGACC      ! accumulation mode

                           ! modal geometric standard deviation:
      REAL( 8 ) :: SGATK      ! Atken mode
      REAL( 8 ) :: SGACC      ! accumulation mode

                           ! natural log of modal geometric standard deviation:
      REAL( 8 ) :: XXLSGAT    ! Aitken mode
      REAL( 8 ) :: XXLSGAC    ! accumulation mode

                           ! coagulation coefficients
      REAL( 8 ) :: QS11, QN11, QS22, QN22,
     &             QS12, QS21, QN12, QV12 

C *** Local Variables:

      INTEGER IBETA, N1, N2A, N2N ! indices for correction factors

      REAL( 8 ) :: I1FM_AT
      REAL( 8 ) :: I1NC_AT
      REAL( 8 ) :: I1_AT
      
      REAL( 8 ) :: I1FM_AC
      REAL( 8 ) :: I1NC_AC
      REAL( 8 ) :: I1_AC
      
      REAL( 8 ) :: I1FM
      REAL( 8 ) :: I1NC
      REAL( 8 ) :: I1
      
      REAL( 8 ) :: CONST II
      
      REAL( 8 ) :: KNGAT, KNGAC

      REAL( 8 ), PARAMETER :: ONE = 1.0D0, TWO = 2.0D0, HALF = 0.5D0
!     REAL( 8 ), PARAMETER :: A = 2.492D0
      REAL( 8 ), PARAMETER :: A = 1.246D0
      REAL( 8 ), PARAMETER :: TWO3RDS = 2.D0 / 3.D0

      REAL( 8 ) :: SQRTTWO  !  sqrt(2.0)
      REAL( 8 ) :: DLGSQT2  !  1/ln( sqrt( 2.0 ) )

      REAL( 8 ) :: ESAT01         ! Aitken mode exp( log^2( sigmag )/8 )
      REAL( 8 ) :: ESAC01         ! accumulation mode exp( log^2( sigmag )/8 )

      REAL( 8 ) :: ESAT04
      REAL( 8 ) :: ESAC04

      REAL( 8 ) :: ESAT05
      REAL( 8 ) :: ESAC05

      REAL( 8 ) :: ESAT08
      REAL( 8 ) :: ESAC08

      REAL( 8 ) :: ESAT09
      REAL( 8 ) :: ESAC09

      REAL( 8 ) :: ESAT16
      REAL( 8 ) :: ESAC16
      
      REAL( 8 ) :: ESAT20
      REAL( 8 ) :: ESAC20
      
      REAL( 8 ) :: ESAT24
      REAL( 8 ) :: ESAC24
     
      REAL( 8 ) :: ESAT25
      REAL( 8 ) :: ESAC25
      
      REAL( 8 ) :: ESAT36
      REAL( 8 ) :: ESAC36

      REAL( 8 ) :: ESAT49
      
      REAL( 8 ) :: ESAT64
      REAL( 8 ) :: ESAC64
            
      REAL( 8 ) :: ESAT100
          
      REAL( 8 ) :: DGAT2, DGAC2, DGAT3, DGAC3
      REAL( 8 ) :: SQDGAT, SQDGAC
      REAL( 8 ) :: SQDGAT5, SQDGAC5
      REAL( 8 ) :: SQDGAT7
      REAL( 8 ) :: R, R2, R3, R4, R5, R6, R8
      REAL( 8 ) :: RI1, RI2, RI3, RI4
      REAL( 8 ) :: RAT
      REAL( 8 ) :: COAGFM0, COAGNC0
      REAL( 8 ) :: COAGFM3, COAGNC3
      REAL( 8 ) :: COAGFM_AT, COAGFM_AC
      REAL( 8 ) :: COAGNC_AT, COAGNC_AC
      REAL( 8 ) :: COAGATAT0
      REAL( 8 ) :: COAGACAC0
      REAL( 8 ) :: COAGATAT2
      REAL( 8 ) :: COAGACAC2
      REAL( 8 ) :: COAGATAC0, COAGATAC3
      REAL( 8 ) :: COAGATAC2
      REAL( 8 ) :: COAGACAT2
      REAL( 8 ) :: XM2AT, XM3AT, XM2AC, XM3AC

C *** correction factors for coagulation rates      
      REAL :: BM0    ( 10 )        ! M0 INTRAmodal FM - RPM values
      REAL :: BM0IJ  ( 10,10,10 )  ! M0 INTERmodal FM
      REAL :: BM3I   ( 10,10,10 )  ! M3 INTERmodal FM- RPM values
      REAL :: BM2II  ( 10 )        ! M2 INTRAmodal FM
      REAL :: BM2IITT( 10 )        ! M2 INTRAmodal total
      REAL :: BM2IJ  ( 10,10,10 )  ! M2 INTERmodal FM i to j
      REAL :: BM2JI  ( 10,10,10 )  ! M2 total INTERmodal  j from i

C *** populate the arrays for the correction factors.

C RPM 0th moment correction factors for UNIMODAL FM coagulation  rates
      DATA BM0 /
     &      0.707106785165097, 0.726148960080488, 0.766430744110958,
     &      0.814106389441342, 0.861679526483207, 0.903600509090092,
     &      0.936578814219156, 0.960098926735545, 0.975646823342881,
     &      0.985397173215326 /

C FSB new FM correction factors for M0 INTERmodal coagulation

      DATA ( BM0IJ(  1,  1,IBETA ), IBETA = 1, 10 ) /
     &  0.628539,  0.639610,  0.664514,  0.696278,  0.731558,
     &  0.768211,  0.804480,  0.838830,  0.870024,  0.897248/
      DATA ( BM0IJ(  1,  2,IBETA ), IBETA = 1, 10 ) /
     &  0.639178,  0.649966,  0.674432,  0.705794,  0.740642,
     &  0.776751,  0.812323,  0.845827,  0.876076,  0.902324/
      DATA ( BM0IJ(  1,  3,IBETA ), IBETA = 1, 10 ) /
     &  0.663109,  0.673464,  0.697147,  0.727637,  0.761425,
     &  0.796155,  0.829978,  0.861419,  0.889424,  0.913417/
      DATA ( BM0IJ(  1,  4,IBETA ), IBETA = 1, 10 ) /
     &  0.693693,  0.703654,  0.726478,  0.755786,  0.787980,
     &  0.820626,  0.851898,  0.880459,  0.905465,  0.926552/
      DATA ( BM0IJ(  1,  5,IBETA ), IBETA = 1, 10 ) /
     &  0.727803,  0.737349,  0.759140,  0.786870,  0.816901,
     &  0.846813,  0.874906,  0.900060,  0.921679,  0.939614/
      DATA ( BM0IJ(  1,  6,IBETA ), IBETA = 1, 10 ) /
     &  0.763461,  0.772483,  0.792930,  0.818599,  0.845905,
     &  0.872550,  0.897051,  0.918552,  0.936701,  0.951528/
      DATA ( BM0IJ(  1,  7,IBETA ), IBETA = 1, 10 ) /
     &  0.799021,  0.807365,  0.826094,  0.849230,  0.873358,
     &  0.896406,  0.917161,  0.935031,  0.949868,  0.961828/
      DATA ( BM0IJ(  1,  8,IBETA ), IBETA = 1, 10 ) /
     &  0.833004,  0.840514,  0.857192,  0.877446,  0.898147,
     &  0.917518,  0.934627,  0.949106,  0.960958,  0.970403/
      DATA ( BM0IJ(  1,  9,IBETA ), IBETA = 1, 10 ) /
     &  0.864172,  0.870734,  0.885153,  0.902373,  0.919640,
     &  0.935494,  0.949257,  0.960733,  0.970016,  0.977346/
      DATA ( BM0IJ(  1, 10,IBETA ), IBETA = 1, 10 ) /
     &  0.891658,  0.897227,  0.909343,  0.923588,  0.937629,
     &  0.950307,  0.961151,  0.970082,  0.977236,  0.982844/
      DATA ( BM0IJ(  2,  1,IBETA ), IBETA = 1, 10 ) /
     &  0.658724,  0.670587,  0.697539,  0.731890,  0.769467,
     &  0.807391,  0.843410,  0.875847,  0.903700,  0.926645/
      DATA ( BM0IJ(  2,  2,IBETA ), IBETA = 1, 10 ) /
     &  0.667070,  0.678820,  0.705538,  0.739591,  0.776758,
     &  0.814118,  0.849415,  0.881020,  0.908006,  0.930121/
      DATA ( BM0IJ(  2,  3,IBETA ), IBETA = 1, 10 ) /
     &  0.686356,  0.697839,  0.723997,  0.757285,  0.793389,
     &  0.829313,  0.862835,  0.892459,  0.917432,  0.937663/
      DATA ( BM0IJ(  2,  4,IBETA ), IBETA = 1, 10 ) /
     &  0.711425,  0.722572,  0.747941,  0.780055,  0.814518,
     &  0.848315,  0.879335,  0.906290,  0.928658,  0.946526/
      DATA ( BM0IJ(  2,  5,IBETA ), IBETA = 1, 10 ) /
     &  0.739575,  0.750307,  0.774633,  0.805138,  0.837408,
     &  0.868504,  0.896517,  0.920421,  0.939932,  0.955299/
      DATA ( BM0IJ(  2,  6,IBETA ), IBETA = 1, 10 ) /
     &  0.769143,  0.779346,  0.802314,  0.830752,  0.860333,
     &  0.888300,  0.913014,  0.933727,  0.950370,  0.963306/
      DATA ( BM0IJ(  2,  7,IBETA ), IBETA = 1, 10 ) /
     &  0.798900,  0.808431,  0.829700,  0.855653,  0.882163,
     &  0.906749,  0.928075,  0.945654,  0.959579,  0.970280/
      DATA ( BM0IJ(  2,  8,IBETA ), IBETA = 1, 10 ) /
     &  0.827826,  0.836542,  0.855808,  0.878954,  0.902174,
     &  0.923316,  0.941345,  0.955989,  0.967450,  0.976174/
      DATA ( BM0IJ(  2,  9,IBETA ), IBETA = 1, 10 ) /
     &  0.855068,  0.862856,  0.879900,  0.900068,  0.919956,
     &  0.937764,  0.952725,  0.964726,  0.974027,  0.981053/
      DATA ( BM0IJ(  2, 10,IBETA ), IBETA = 1, 10 ) /
     &  0.879961,  0.886755,  0.901484,  0.918665,  0.935346,
     &  0.950065,  0.962277,  0.971974,  0.979432,  0.985033/
      DATA ( BM0IJ(  3,  1,IBETA ), IBETA = 1, 10 ) /
     &  0.724166,  0.735474,  0.761359,  0.794045,  0.828702,
     &  0.862061,  0.891995,  0.917385,  0.937959,  0.954036/
      DATA ( BM0IJ(  3,  2,IBETA ), IBETA = 1, 10 ) /
     &  0.730416,  0.741780,  0.767647,  0.800116,  0.834344,
     &  0.867093,  0.896302,  0.920934,  0.940790,  0.956237/
      DATA ( BM0IJ(  3,  3,IBETA ), IBETA = 1, 10 ) /
     &  0.745327,  0.756664,  0.782255,  0.814026,  0.847107,
     &  0.878339,  0.905820,  0.928699,  0.946931,  0.960977/
      DATA ( BM0IJ(  3,  4,IBETA ), IBETA = 1, 10 ) /
     &  0.765195,  0.776312,  0.801216,  0.831758,  0.863079,
     &  0.892159,  0.917319,  0.937939,  0.954145,  0.966486/
      DATA ( BM0IJ(  3,  5,IBETA ), IBETA = 1, 10 ) /
     &  0.787632,  0.798347,  0.822165,  0.850985,  0.880049,
     &  0.906544,  0.929062,  0.947218,  0.961288,  0.971878/
      DATA ( BM0IJ(  3,  6,IBETA ), IBETA = 1, 10 ) /
     &  0.811024,  0.821179,  0.843557,  0.870247,  0.896694,
     &  0.920365,  0.940131,  0.955821,  0.967820,  0.976753/
      DATA ( BM0IJ(  3,  7,IBETA ), IBETA = 1, 10 ) /
     &  0.834254,  0.843709,  0.864356,  0.888619,  0.912245,
     &  0.933019,  0.950084,  0.963438,  0.973530,  0.980973/
      DATA ( BM0IJ(  3,  8,IBETA ), IBETA = 1, 10 ) /
     &  0.856531,  0.865176,  0.883881,  0.905544,  0.926290,
     &  0.944236,  0.958762,  0.969988,  0.978386,  0.984530/
      DATA ( BM0IJ(  3,  9,IBETA ), IBETA = 1, 10 ) /
     &  0.877307,  0.885070,  0.901716,  0.920729,  0.938663,
     &  0.953951,  0.966169,  0.975512,  0.982442,  0.987477/
      DATA ( BM0IJ(  3, 10,IBETA ), IBETA = 1, 10 ) /
     &  0.896234,  0.903082,  0.917645,  0.934069,  0.949354,
     &  0.962222,  0.972396,  0.980107,  0.985788,  0.989894/
      DATA ( BM0IJ(  4,  1,IBETA ), IBETA = 1, 10 ) /
     &  0.799294,  0.809144,  0.831293,  0.858395,  0.885897,
     &  0.911031,  0.932406,  0.949642,  0.963001,  0.973062/
      DATA ( BM0IJ(  4,  2,IBETA ), IBETA = 1, 10 ) /
     &  0.804239,  0.814102,  0.836169,  0.862984,  0.890003,
     &  0.914535,  0.935274,  0.951910,  0.964748,  0.974381/
      DATA ( BM0IJ(  4,  3,IBETA ), IBETA = 1, 10 ) /
     &  0.815910,  0.825708,  0.847403,  0.873389,  0.899185,
     &  0.922275,  0.941543,  0.956826,  0.968507,  0.977204/
      DATA ( BM0IJ(  4,  4,IBETA ), IBETA = 1, 10 ) /
     &  0.831348,  0.840892,  0.861793,  0.886428,  0.910463,
     &  0.931614,  0.948993,  0.962593,  0.972872,  0.980456/
      DATA ( BM0IJ(  4,  5,IBETA ), IBETA = 1, 10 ) /
     &  0.848597,  0.857693,  0.877402,  0.900265,  0.922180,
     &  0.941134,  0.956464,  0.968298,  0.977143,  0.983611/
      DATA ( BM0IJ(  4,  6,IBETA ), IBETA = 1, 10 ) /
     &  0.866271,  0.874764,  0.892984,  0.913796,  0.933407,
     &  0.950088,  0.963380,  0.973512,  0.981006,  0.986440/
      DATA ( BM0IJ(  4,  7,IBETA ), IBETA = 1, 10 ) /
     &  0.883430,  0.891216,  0.907762,  0.926388,  0.943660,
     &  0.958127,  0.969499,  0.978070,  0.984351,  0.988872/
      DATA ( BM0IJ(  4,  8,IBETA ), IBETA = 1, 10 ) /
     &  0.899483,  0.906505,  0.921294,  0.937719,  0.952729,
     &  0.965131,  0.974762,  0.981950,  0.987175,  0.990912/
      DATA ( BM0IJ(  4,  9,IBETA ), IBETA = 1, 10 ) /
     &  0.914096,  0.920337,  0.933373,  0.947677,  0.960579,
     &  0.971111,  0.979206,  0.985196,  0.989520,  0.992597/
      DATA ( BM0IJ(  4, 10,IBETA ), IBETA = 1, 10 ) /
     &  0.927122,  0.932597,  0.943952,  0.956277,  0.967268,
     &  0.976147,  0.982912,  0.987882,  0.991450,  0.993976/
      DATA ( BM0IJ(  5,  1,IBETA ), IBETA = 1, 10 ) /
     &  0.865049,  0.872851,  0.889900,  0.909907,  0.929290,
     &  0.946205,  0.959991,  0.970706,  0.978764,  0.984692/
      DATA ( BM0IJ(  5,  2,IBETA ), IBETA = 1, 10 ) /
     &  0.868989,  0.876713,  0.893538,  0.913173,  0.932080,
     &  0.948484,  0.961785,  0.972080,  0.979796,  0.985457/
      DATA ( BM0IJ(  5,  3,IBETA ), IBETA = 1, 10 ) /
     &  0.878010,  0.885524,  0.901756,  0.920464,  0.938235,
     &  0.953461,  0.965672,  0.975037,  0.982005,  0.987085/
      DATA ( BM0IJ(  5,  4,IBETA ), IBETA = 1, 10 ) /
     &  0.889534,  0.896698,  0.912012,  0.929395,  0.945647,
     &  0.959366,  0.970227,  0.978469,  0.984547,  0.988950/
      DATA ( BM0IJ(  5,  5,IBETA ), IBETA = 1, 10 ) /
     &  0.902033,  0.908713,  0.922848,  0.938648,  0.953186,
     &  0.965278,  0.974729,  0.981824,  0.987013,  0.990746/
      DATA ( BM0IJ(  5,  6,IBETA ), IBETA = 1, 10 ) /
     &  0.914496,  0.920599,  0.933389,  0.947485,  0.960262,
     &  0.970743,  0.978839,  0.984858,  0.989225,  0.992348/
      DATA ( BM0IJ(  5,  7,IBETA ), IBETA = 1, 10 ) /
     &  0.926281,  0.931761,  0.943142,  0.955526,  0.966600,
     &  0.975573,  0.982431,  0.987485,  0.991128,  0.993718/
      DATA ( BM0IJ(  5,  8,IBETA ), IBETA = 1, 10 ) /
     &  0.937029,  0.941877,  0.951868,  0.962615,  0.972112,
     &  0.979723,  0.985488,  0.989705,  0.992725,  0.994863/
      DATA ( BM0IJ(  5,  9,IBETA ), IBETA = 1, 10 ) /
     &  0.946580,  0.950819,  0.959494,  0.968732,  0.976811,
     &  0.983226,  0.988047,  0.991550,  0.994047,  0.995806/
      DATA ( BM0IJ(  5, 10,IBETA ), IBETA = 1, 10 ) /
     &  0.954909,  0.958581,  0.966049,  0.973933,  0.980766,
     &  0.986149,  0.990166,  0.993070,  0.995130,  0.996577/
      DATA ( BM0IJ(  6,  1,IBETA ), IBETA = 1, 10 ) /
     &  0.914182,  0.919824,  0.931832,  0.945387,  0.957999,
     &  0.968606,  0.976982,  0.983331,  0.988013,  0.991407/
      DATA ( BM0IJ(  6,  2,IBETA ), IBETA = 1, 10 ) /
     &  0.917139,  0.922665,  0.934395,  0.947580,  0.959792,
     &  0.970017,  0.978062,  0.984138,  0.988609,  0.991843/
      DATA ( BM0IJ(  6,  3,IBETA ), IBETA = 1, 10 ) /
     &  0.923742,  0.928990,  0.940064,  0.952396,  0.963699,
     &  0.973070,  0.980381,  0.985866,  0.989878,  0.992768/
      DATA ( BM0IJ(  6,  4,IBETA ), IBETA = 1, 10 ) /
     &  0.931870,  0.936743,  0.946941,  0.958162,  0.968318,
     &  0.976640,  0.983069,  0.987853,  0.991330,  0.993822/
      DATA ( BM0IJ(  6,  5,IBETA ), IBETA = 1, 10 ) /
     &  0.940376,  0.944807,  0.954004,  0.963999,  0.972928,
     &  0.980162,  0.985695,  0.989779,  0.992729,  0.994833/
      DATA ( BM0IJ(  6,  6,IBETA ), IBETA = 1, 10 ) /
     &  0.948597,  0.952555,  0.960703,  0.969454,  0.977181,
     &  0.983373,  0.988067,  0.991507,  0.993977,  0.995730/
      DATA ( BM0IJ(  6,  7,IBETA ), IBETA = 1, 10 ) /
     &  0.956167,  0.959648,  0.966763,  0.974326,  0.980933,
     &  0.986177,  0.990121,  0.992993,  0.995045,  0.996495/
      DATA ( BM0IJ(  6,  8,IBETA ), IBETA = 1, 10 ) /
     &  0.962913,  0.965937,  0.972080,  0.978552,  0.984153,
     &  0.988563,  0.991857,  0.994242,  0.995938,  0.997133/
      DATA ( BM0IJ(  6,  9,IBETA ), IBETA = 1, 10 ) /
     &  0.968787,  0.971391,  0.976651,  0.982148,  0.986869,
     &  0.990560,  0.993301,  0.995275,  0.996675,  0.997657/
      DATA ( BM0IJ(  6, 10,IBETA ), IBETA = 1, 10 ) /
     &  0.973822,  0.976047,  0.980523,  0.985170,  0.989134,
     &  0.992215,  0.994491,  0.996124,  0.997277,  0.998085/
      DATA ( BM0IJ(  7,  1,IBETA ), IBETA = 1, 10 ) /
     &  0.947410,  0.951207,  0.959119,  0.967781,  0.975592,
     &  0.981981,  0.986915,  0.990590,  0.993266,  0.995187/
      DATA ( BM0IJ(  7,  2,IBETA ), IBETA = 1, 10 ) /
     &  0.949477,  0.953161,  0.960824,  0.969187,  0.976702,
     &  0.982831,  0.987550,  0.991057,  0.993606,  0.995434/
      DATA ( BM0IJ(  7,  3,IBETA ), IBETA = 1, 10 ) /
     &  0.954008,  0.957438,  0.964537,  0.972232,  0.979095,
     &  0.984653,  0.988907,  0.992053,  0.994330,  0.995958/
      DATA ( BM0IJ(  7,  4,IBETA ), IBETA = 1, 10 ) /
     &  0.959431,  0.962539,  0.968935,  0.975808,  0.981882,
     &  0.986759,  0.990466,  0.993190,  0.995153,  0.996552/
      DATA ( BM0IJ(  7,  5,IBETA ), IBETA = 1, 10 ) /
     &  0.964932,  0.967693,  0.973342,  0.979355,  0.984620,
     &  0.988812,  0.991974,  0.994285,  0.995943,  0.997119/
      DATA ( BM0IJ(  7,  6,IBETA ), IBETA = 1, 10 ) /
     &  0.970101,  0.972517,  0.977428,  0.982612,  0.987110,
     &  0.990663,  0.993326,  0.995261,  0.996644,  0.997621/
      DATA ( BM0IJ(  7,  7,IBETA ), IBETA = 1, 10 ) /
     &  0.974746,  0.976834,  0.981055,  0.985475,  0.989280,
     &  0.992265,  0.994488,  0.996097,  0.997241,  0.998048/
      DATA ( BM0IJ(  7,  8,IBETA ), IBETA = 1, 10 ) /
     &  0.978804,  0.980591,  0.984187,  0.987927,  0.991124,
     &  0.993617,  0.995464,  0.996795,  0.997739,  0.998403/
      DATA ( BM0IJ(  7,  9,IBETA ), IBETA = 1, 10 ) /
     &  0.982280,  0.983799,  0.986844,  0.989991,  0.992667,
     &  0.994742,  0.996273,  0.997372,  0.998149,  0.998695/
      DATA ( BM0IJ(  7, 10,IBETA ), IBETA = 1, 10 ) /
     &  0.985218,  0.986503,  0.989071,  0.991711,  0.993945,
     &  0.995669,  0.996937,  0.997844,  0.998484,  0.998932/
      DATA ( BM0IJ(  8,  1,IBETA ), IBETA = 1, 10 ) /
     &  0.968507,  0.970935,  0.975916,  0.981248,  0.985947,
     &  0.989716,  0.992580,  0.994689,  0.996210,  0.997297/
      DATA ( BM0IJ(  8,  2,IBETA ), IBETA = 1, 10 ) /
     &  0.969870,  0.972210,  0.977002,  0.982119,  0.986619,
     &  0.990219,  0.992951,  0.994958,  0.996405,  0.997437/
      DATA ( BM0IJ(  8,  3,IBETA ), IBETA = 1, 10 ) /
     &  0.972820,  0.974963,  0.979339,  0.983988,  0.988054,
     &  0.991292,  0.993738,  0.995529,  0.996817,  0.997734/
      DATA ( BM0IJ(  8,  4,IBETA ), IBETA = 1, 10 ) /
     &  0.976280,  0.978186,  0.982060,  0.986151,  0.989706,
     &  0.992520,  0.994636,  0.996179,  0.997284,  0.998069/
      DATA ( BM0IJ(  8,  5,IBETA ), IBETA = 1, 10 ) /
     &  0.979711,  0.981372,  0.984735,  0.988263,  0.991309,
     &  0.993706,  0.995499,  0.996801,  0.997730,  0.998389/
      DATA ( BM0IJ(  8,  6,IBETA ), IBETA = 1, 10 ) /
     &  0.982863,  0.984292,  0.987172,  0.990174,  0.992750,
     &  0.994766,  0.996266,  0.997352,  0.998125,  0.998670/
      DATA ( BM0IJ(  8,  7,IBETA ), IBETA = 1, 10 ) /
     &  0.985642,  0.986858,  0.989301,  0.991834,  0.993994,
     &  0.995676,  0.996923,  0.997822,  0.998460,  0.998910/
      DATA ( BM0IJ(  8,  8,IBETA ), IBETA = 1, 10 ) /
     &  0.988029,  0.989058,  0.991116,  0.993240,  0.995043,
     &  0.996440,  0.997472,  0.998214,  0.998739,  0.999108/
      DATA ( BM0IJ(  8,  9,IBETA ), IBETA = 1, 10 ) /
     &  0.990046,  0.990912,  0.992640,  0.994415,  0.995914,
     &  0.997073,  0.997925,  0.998536,  0.998968,  0.999271/
      DATA ( BM0IJ(  8, 10,IBETA ), IBETA = 1, 10 ) /
     &  0.991732,  0.992459,  0.993906,  0.995386,  0.996633,
     &  0.997592,  0.998296,  0.998799,  0.999154,  0.999403/
      DATA ( BM0IJ(  9,  1,IBETA ), IBETA = 1, 10 ) /
     &  0.981392,  0.982893,  0.985938,  0.989146,  0.991928,
     &  0.994129,  0.995783,  0.996991,  0.997857,  0.998473/
      DATA ( BM0IJ(  9,  2,IBETA ), IBETA = 1, 10 ) /
     &  0.982254,  0.983693,  0.986608,  0.989673,  0.992328,
     &  0.994424,  0.995998,  0.997146,  0.997969,  0.998553/
      DATA ( BM0IJ(  9,  3,IBETA ), IBETA = 1, 10 ) /
     &  0.984104,  0.985407,  0.988040,  0.990798,  0.993178,
     &  0.995052,  0.996454,  0.997474,  0.998204,  0.998722/
      DATA ( BM0IJ(  9,  4,IBETA ), IBETA = 1, 10 ) /
     &  0.986243,  0.987386,  0.989687,  0.992087,  0.994149,
     &  0.995765,  0.996971,  0.997846,  0.998470,  0.998913/
      DATA ( BM0IJ(  9,  5,IBETA ), IBETA = 1, 10 ) /
     &  0.988332,  0.989313,  0.991284,  0.993332,  0.995082,
     &  0.996449,  0.997465,  0.998200,  0.998723,  0.999093/
      DATA ( BM0IJ(  9,  6,IBETA ), IBETA = 1, 10 ) /
     &  0.990220,  0.991053,  0.992721,  0.994445,  0.995914,
     &  0.997056,  0.997902,  0.998513,  0.998947,  0.999253/
      DATA ( BM0IJ(  9,  7,IBETA ), IBETA = 1, 10 ) /
     &  0.991859,  0.992561,  0.993961,  0.995403,  0.996626,
     &  0.997574,  0.998274,  0.998778,  0.999136,  0.999387/
      DATA ( BM0IJ(  9,  8,IBETA ), IBETA = 1, 10 ) /
     &  0.993250,  0.993837,  0.995007,  0.996208,  0.997223,
     &  0.998007,  0.998584,  0.998999,  0.999293,  0.999499/
      DATA ( BM0IJ(  9,  9,IBETA ), IBETA = 1, 10 ) /
     &  0.994413,  0.994903,  0.995878,  0.996876,  0.997716,
     &  0.998363,  0.998839,  0.999180,  0.999421,  0.999591/
      DATA ( BM0IJ(  9, 10,IBETA ), IBETA = 1, 10 ) /
     &  0.995376,  0.995785,  0.996597,  0.997425,  0.998121,
     &  0.998655,  0.999048,  0.999328,  0.999526,  0.999665/
      DATA ( BM0IJ( 10,  1,IBETA ), IBETA = 1, 10 ) /
     &  0.989082,  0.989991,  0.991819,  0.993723,  0.995357,
     &  0.996637,  0.997592,  0.998286,  0.998781,  0.999132/
      DATA ( BM0IJ( 10,  2,IBETA ), IBETA = 1, 10 ) /
     &  0.989613,  0.990480,  0.992224,  0.994039,  0.995594,
     &  0.996810,  0.997717,  0.998375,  0.998845,  0.999178/
      DATA ( BM0IJ( 10,  3,IBETA ), IBETA = 1, 10 ) /
     &  0.990744,  0.991523,  0.993086,  0.994708,  0.996094,
     &  0.997176,  0.997981,  0.998564,  0.998980,  0.999274/
      DATA ( BM0IJ( 10,  4,IBETA ), IBETA = 1, 10 ) /
     &  0.992041,  0.992716,  0.994070,  0.995470,  0.996662,
     &  0.997591,  0.998280,  0.998778,  0.999133,  0.999383/
      DATA ( BM0IJ( 10,  5,IBETA ), IBETA = 1, 10 ) /
     &  0.993292,  0.993867,  0.995015,  0.996199,  0.997205,
     &  0.997985,  0.998564,  0.998981,  0.999277,  0.999487/
      DATA ( BM0IJ( 10,  6,IBETA ), IBETA = 1, 10 ) /
     &  0.994411,  0.994894,  0.995857,  0.996847,  0.997685,
     &  0.998334,  0.998814,  0.999159,  0.999404,  0.999577/
      DATA ( BM0IJ( 10,  7,IBETA ), IBETA = 1, 10 ) /
     &  0.995373,  0.995776,  0.996577,  0.997400,  0.998094,
     &  0.998630,  0.999026,  0.999310,  0.999512,  0.999654/
      DATA ( BM0IJ( 10,  8,IBETA ), IBETA = 1, 10 ) /
     &  0.996181,  0.996516,  0.997181,  0.997861,  0.998435,
     &  0.998877,  0.999202,  0.999435,  0.999601,  0.999717/
      DATA ( BM0IJ( 10,  9,IBETA ), IBETA = 1, 10 ) /
     &  0.996851,  0.997128,  0.997680,  0.998242,  0.998715,
     &  0.999079,  0.999346,  0.999538,  0.999673,  0.999769/
      DATA ( BM0IJ( 10, 10,IBETA ), IBETA = 1, 10 ) /
     &  0.997402,  0.997632,  0.998089,  0.998554,  0.998945,
     &  0.999244,  0.999464,  0.999622,  0.999733,  0.999811/

C RPM....   3rd moment nuclei mode corr. fac. for bimodal FM coag rate
       DATA ( BM3I(  1, 1,IBETA ), IBETA = 1, 10 ) /
     + 0.70708,0.71681,0.73821,0.76477,0.79350,0.82265,0.85090,0.87717,
     + 0.90069,0.92097/
       DATA ( BM3I(  1, 2,IBETA ), IBETA = 1, 10 ) /
     + 0.72172,0.73022,0.74927,0.77324,0.79936,0.82601,0.85199,0.87637,
     + 0.89843,0.91774/
       DATA ( BM3I(  1, 3,IBETA ), IBETA = 1, 10 ) /
     + 0.78291,0.78896,0.80286,0.82070,0.84022,0.85997,0.87901,0.89669,
     + 0.91258,0.92647/
       DATA ( BM3I(  1, 4,IBETA ), IBETA = 1, 10 ) /
     + 0.87760,0.88147,0.89025,0.90127,0.91291,0.92420,0.93452,0.94355,
     + 0.95113,0.95726/
       DATA ( BM3I(  1, 5,IBETA ), IBETA = 1, 10 ) /
     + 0.94988,0.95184,0.95612,0.96122,0.96628,0.97085,0.97467,0.97763,
     + 0.97971,0.98089/
       DATA ( BM3I(  1, 6,IBETA ), IBETA = 1, 10 ) /
     + 0.98318,0.98393,0.98551,0.98728,0.98889,0.99014,0.99095,0.99124,
     + 0.99100,0.99020/
       DATA ( BM3I(  1, 7,IBETA ), IBETA = 1, 10 ) /
     + 0.99480,0.99504,0.99551,0.99598,0.99629,0.99635,0.99611,0.99550,
     + 0.99450,0.99306/
       DATA ( BM3I(  1, 8,IBETA ), IBETA = 1, 10 ) /
     + 0.99842,0.99848,0.99858,0.99861,0.99850,0.99819,0.99762,0.99674,
     + 0.99550,0.99388/
       DATA ( BM3I(  1, 9,IBETA ), IBETA = 1, 10 ) /
     + 0.99951,0.99951,0.99949,0.99939,0.99915,0.99872,0.99805,0.99709,
     + 0.99579,0.99411/
       DATA ( BM3I(  1,10,IBETA ), IBETA = 1, 10 ) /
     + 0.99984,0.99982,0.99976,0.99962,0.99934,0.99888,0.99818,0.99719,
     + 0.99587,0.99417/
       DATA ( BM3I(  2, 1,IBETA ), IBETA = 1, 10 ) /
     + 0.72957,0.73993,0.76303,0.79178,0.82245,0.85270,0.88085,0.90578,
     + 0.92691,0.94415/
       DATA ( BM3I(  2, 2,IBETA ), IBETA = 1, 10 ) /
     + 0.72319,0.73320,0.75547,0.78323,0.81307,0.84287,0.87107,0.89651,
     + 0.91852,0.93683/
       DATA ( BM3I(  2, 3,IBETA ), IBETA = 1, 10 ) /
     + 0.74413,0.75205,0.76998,0.79269,0.81746,0.84258,0.86685,0.88938,
     + 0.90953,0.92695/
       DATA ( BM3I(  2, 4,IBETA ), IBETA = 1, 10 ) /
     + 0.82588,0.83113,0.84309,0.85825,0.87456,0.89072,0.90594,0.91972,
     + 0.93178,0.94203/
       DATA ( BM3I(  2, 5,IBETA ), IBETA = 1, 10 ) /
     + 0.91886,0.92179,0.92831,0.93624,0.94434,0.95192,0.95856,0.96409,
     + 0.96845,0.97164/
       DATA ( BM3I(  2, 6,IBETA ), IBETA = 1, 10 ) /
     + 0.97129,0.97252,0.97515,0.97818,0.98108,0.98354,0.98542,0.98665,
     + 0.98721,0.98709/
       DATA ( BM3I(  2, 7,IBETA ), IBETA = 1, 10 ) /
     + 0.99104,0.99145,0.99230,0.99320,0.99394,0.99439,0.99448,0.99416,
     + 0.99340,0.99217/
       DATA ( BM3I(  2, 8,IBETA ), IBETA = 1, 10 ) /
     + 0.99730,0.99741,0.99763,0.99779,0.99782,0.99762,0.99715,0.99636,
     + 0.99519,0.99363/
       DATA ( BM3I(  2, 9,IBETA ), IBETA = 1, 10 ) /
     + 0.99917,0.99919,0.99921,0.99915,0.99895,0.99856,0.99792,0.99698,
     + 0.99570,0.99404/
       DATA ( BM3I(  2,10,IBETA ), IBETA = 1, 10 ) /
     + 0.99973,0.99973,0.99968,0.99955,0.99928,0.99883,0.99814,0.99716,
     + 0.99584,0.99415/
       DATA ( BM3I(  3, 1,IBETA ), IBETA = 1, 10 ) /
     + 0.78358,0.79304,0.81445,0.84105,0.86873,0.89491,0.91805,0.93743,
     + 0.95300,0.96510/
       DATA ( BM3I(  3, 2,IBETA ), IBETA = 1, 10 ) /
     + 0.76412,0.77404,0.79635,0.82404,0.85312,0.88101,0.90610,0.92751,
     + 0.94500,0.95879/
       DATA ( BM3I(  3, 3,IBETA ), IBETA = 1, 10 ) /
     + 0.74239,0.75182,0.77301,0.79956,0.82809,0.85639,0.88291,0.90658,
     + 0.92683,0.94350/
       DATA ( BM3I(  3, 4,IBETA ), IBETA = 1, 10 ) /
     + 0.78072,0.78758,0.80317,0.82293,0.84437,0.86589,0.88643,0.90526,
     + 0.92194,0.93625/
       DATA ( BM3I(  3, 5,IBETA ), IBETA = 1, 10 ) /
     + 0.87627,0.88044,0.88981,0.90142,0.91357,0.92524,0.93585,0.94510,
     + 0.95285,0.95911/
       DATA ( BM3I(  3, 6,IBETA ), IBETA = 1, 10 ) /
     + 0.95176,0.95371,0.95796,0.96297,0.96792,0.97233,0.97599,0.97880,
     + 0.98072,0.98178/
       DATA ( BM3I(  3, 7,IBETA ), IBETA = 1, 10 ) /
     + 0.98453,0.98523,0.98670,0.98833,0.98980,0.99092,0.99160,0.99179,
     + 0.99145,0.99058/
       DATA ( BM3I(  3, 8,IBETA ), IBETA = 1, 10 ) /
     + 0.99534,0.99555,0.99597,0.99637,0.99662,0.99663,0.99633,0.99569,
     + 0.99465,0.99318/
       DATA ( BM3I(  3, 9,IBETA ), IBETA = 1, 10 ) /
     + 0.99859,0.99864,0.99872,0.99873,0.99860,0.99827,0.99768,0.99679,
     + 0.99555,0.99391/
       DATA ( BM3I(  3,10,IBETA ), IBETA = 1, 10 ) /
     + 0.99956,0.99956,0.99953,0.99942,0.99918,0.99875,0.99807,0.99711,
     + 0.99580,0.99412/
       DATA ( BM3I(  4, 1,IBETA ), IBETA = 1, 10 ) /
     + 0.84432,0.85223,0.86990,0.89131,0.91280,0.93223,0.94861,0.96172,
     + 0.97185,0.97945/
       DATA ( BM3I(  4, 2,IBETA ), IBETA = 1, 10 ) /
     + 0.82299,0.83164,0.85101,0.87463,0.89857,0.92050,0.93923,0.95443,
     + 0.96629,0.97529/
       DATA ( BM3I(  4, 3,IBETA ), IBETA = 1, 10 ) /
     + 0.77870,0.78840,0.81011,0.83690,0.86477,0.89124,0.91476,0.93460,
     + 0.95063,0.96316/
       DATA ( BM3I(  4, 4,IBETA ), IBETA = 1, 10 ) /
     + 0.76386,0.77233,0.79147,0.81557,0.84149,0.86719,0.89126,0.91275,
     + 0.93116,0.94637/
       DATA ( BM3I(  4, 5,IBETA ), IBETA = 1, 10 ) /
     + 0.82927,0.83488,0.84756,0.86346,0.88040,0.89704,0.91257,0.92649,
     + 0.93857,0.94874/
       DATA ( BM3I(  4, 6,IBETA ), IBETA = 1, 10 ) /
     + 0.92184,0.92481,0.93136,0.93925,0.94724,0.95462,0.96104,0.96634,
     + 0.97048,0.97348/
       DATA ( BM3I(  4, 7,IBETA ), IBETA = 1, 10 ) /
     + 0.97341,0.97457,0.97706,0.97991,0.98260,0.98485,0.98654,0.98760,
     + 0.98801,0.98777/
       DATA ( BM3I(  4, 8,IBETA ), IBETA = 1, 10 ) /
     + 0.99192,0.99229,0.99305,0.99385,0.99449,0.99486,0.99487,0.99449,
     + 0.99367,0.99239/
       DATA ( BM3I(  4, 9,IBETA ), IBETA = 1, 10 ) /
     + 0.99758,0.99768,0.99787,0.99800,0.99799,0.99777,0.99727,0.99645,
     + 0.99527,0.99369/
       DATA ( BM3I(  4,10,IBETA ), IBETA = 1, 10 ) /
     + 0.99926,0.99928,0.99928,0.99921,0.99900,0.99860,0.99795,0.99701,
     + 0.99572,0.99405/
       DATA ( BM3I(  5, 1,IBETA ), IBETA = 1, 10 ) /
     + 0.89577,0.90190,0.91522,0.93076,0.94575,0.95876,0.96932,0.97751,
     + 0.98367,0.98820/
       DATA ( BM3I(  5, 2,IBETA ), IBETA = 1, 10 ) /
     + 0.87860,0.88547,0.90052,0.91828,0.93557,0.95075,0.96319,0.97292,
     + 0.98028,0.98572/
       DATA ( BM3I(  5, 3,IBETA ), IBETA = 1, 10 ) /
     + 0.83381,0.84240,0.86141,0.88425,0.90707,0.92770,0.94510,0.95906,
     + 0.96986,0.97798/
       DATA ( BM3I(  5, 4,IBETA ), IBETA = 1, 10 ) /
     + 0.78530,0.79463,0.81550,0.84127,0.86813,0.89367,0.91642,0.93566,
     + 0.95125,0.96347/
       DATA ( BM3I(  5, 5,IBETA ), IBETA = 1, 10 ) /
     + 0.79614,0.80332,0.81957,0.84001,0.86190,0.88351,0.90368,0.92169,
     + 0.93718,0.95006/
       DATA ( BM3I(  5, 6,IBETA ), IBETA = 1, 10 ) /
     + 0.88192,0.88617,0.89565,0.90728,0.91931,0.93076,0.94107,0.94997,
     + 0.95739,0.96333/
       DATA ( BM3I(  5, 7,IBETA ), IBETA = 1, 10 ) /
     + 0.95509,0.95698,0.96105,0.96583,0.97048,0.97460,0.97796,0.98050,
     + 0.98218,0.98304/
       DATA ( BM3I(  5, 8,IBETA ), IBETA = 1, 10 ) /
     + 0.98596,0.98660,0.98794,0.98943,0.99074,0.99172,0.99227,0.99235,
     + 0.99192,0.99096/
       DATA ( BM3I(  5, 9,IBETA ), IBETA = 1, 10 ) /
     + 0.99581,0.99600,0.99637,0.99672,0.99691,0.99687,0.99653,0.99585,
     + 0.99478,0.99329/
       DATA ( BM3I(  5,10,IBETA ), IBETA = 1, 10 ) /
     + 0.99873,0.99878,0.99884,0.99883,0.99869,0.99834,0.99774,0.99684,
     + 0.99558,0.99394/
       DATA ( BM3I(  6, 1,IBETA ), IBETA = 1, 10 ) /
     + 0.93335,0.93777,0.94711,0.95764,0.96741,0.97562,0.98210,0.98701,
     + 0.99064,0.99327/
       DATA ( BM3I(  6, 2,IBETA ), IBETA = 1, 10 ) /
     + 0.92142,0.92646,0.93723,0.94947,0.96096,0.97069,0.97842,0.98431,
     + 0.98868,0.99186/
       DATA ( BM3I(  6, 3,IBETA ), IBETA = 1, 10 ) /
     + 0.88678,0.89351,0.90810,0.92508,0.94138,0.95549,0.96693,0.97578,
     + 0.98243,0.98731/
       DATA ( BM3I(  6, 4,IBETA ), IBETA = 1, 10 ) /
     + 0.83249,0.84124,0.86051,0.88357,0.90655,0.92728,0.94477,0.95880,
     + 0.96964,0.97779/
       DATA ( BM3I(  6, 5,IBETA ), IBETA = 1, 10 ) /
     + 0.79593,0.80444,0.82355,0.84725,0.87211,0.89593,0.91735,0.93566,
     + 0.95066,0.96255/
       DATA ( BM3I(  6, 6,IBETA ), IBETA = 1, 10 ) /
     + 0.84124,0.84695,0.85980,0.87575,0.89256,0.90885,0.92383,0.93704,
     + 0.94830,0.95761/
       DATA ( BM3I(  6, 7,IBETA ), IBETA = 1, 10 ) /
     + 0.92721,0.93011,0.93647,0.94406,0.95166,0.95862,0.96460,0.96949,
     + 0.97326,0.97595/
       DATA ( BM3I(  6, 8,IBETA ), IBETA = 1, 10 ) /
     + 0.97573,0.97681,0.97913,0.98175,0.98421,0.98624,0.98772,0.98860,
     + 0.98885,0.98847/
       DATA ( BM3I(  6, 9,IBETA ), IBETA = 1, 10 ) /
     + 0.99271,0.99304,0.99373,0.99444,0.99499,0.99528,0.99522,0.99477,
     + 0.99390,0.99258/
       DATA ( BM3I(  6,10,IBETA ), IBETA = 1, 10 ) /
     + 0.99782,0.99791,0.99807,0.99817,0.99813,0.99788,0.99737,0.99653,
     + 0.99533,0.99374/
       DATA ( BM3I(  7, 1,IBETA ), IBETA = 1, 10 ) /
     + 0.95858,0.96158,0.96780,0.97460,0.98073,0.98575,0.98963,0.99252,
     + 0.99463,0.99615/
       DATA ( BM3I(  7, 2,IBETA ), IBETA = 1, 10 ) /
     + 0.95091,0.95438,0.96163,0.96962,0.97688,0.98286,0.98751,0.99099,
     + 0.99353,0.99536/
       DATA ( BM3I(  7, 3,IBETA ), IBETA = 1, 10 ) /
     + 0.92751,0.93233,0.94255,0.95406,0.96473,0.97366,0.98070,0.98602,
     + 0.98994,0.99278/
       DATA ( BM3I(  7, 4,IBETA ), IBETA = 1, 10 ) /
     + 0.88371,0.89075,0.90595,0.92351,0.94028,0.95474,0.96642,0.97544,
     + 0.98220,0.98715/
       DATA ( BM3I(  7, 5,IBETA ), IBETA = 1, 10 ) /
     + 0.82880,0.83750,0.85671,0.87980,0.90297,0.92404,0.94195,0.95644,
     + 0.96772,0.97625/
       DATA ( BM3I(  7, 6,IBETA ), IBETA = 1, 10 ) /
     + 0.81933,0.82655,0.84279,0.86295,0.88412,0.90449,0.92295,0.93890,
     + 0.95215,0.96281/
       DATA ( BM3I(  7, 7,IBETA ), IBETA = 1, 10 ) /
     + 0.89099,0.89519,0.90448,0.91577,0.92732,0.93820,0.94789,0.95616,
     + 0.96297,0.96838/
       DATA ( BM3I(  7, 8,IBETA ), IBETA = 1, 10 ) /
     + 0.95886,0.96064,0.96448,0.96894,0.97324,0.97701,0.98004,0.98228,
     + 0.98371,0.98435/
       DATA ( BM3I(  7, 9,IBETA ), IBETA = 1, 10 ) /
     + 0.98727,0.98786,0.98908,0.99043,0.99160,0.99245,0.99288,0.99285,
     + 0.99234,0.99131/
       DATA ( BM3I(  7,10,IBETA ), IBETA = 1, 10 ) /
     + 0.99621,0.99638,0.99671,0.99700,0.99715,0.99707,0.99670,0.99599,
     + 0.99489,0.99338/
       DATA ( BM3I(  8, 1,IBETA ), IBETA = 1, 10 ) /
     + 0.97470,0.97666,0.98064,0.98491,0.98867,0.99169,0.99399,0.99569,
     + 0.99691,0.99779/
       DATA ( BM3I(  8, 2,IBETA ), IBETA = 1, 10 ) /
     + 0.96996,0.97225,0.97693,0.98196,0.98643,0.99003,0.99279,0.99482,
     + 0.99630,0.99735/
       DATA ( BM3I(  8, 3,IBETA ), IBETA = 1, 10 ) /
     + 0.95523,0.95848,0.96522,0.97260,0.97925,0.98468,0.98888,0.99200,
     + 0.99427,0.99590/
       DATA ( BM3I(  8, 4,IBETA ), IBETA = 1, 10 ) /
     + 0.92524,0.93030,0.94098,0.95294,0.96397,0.97317,0.98038,0.98582,
     + 0.98981,0.99270/
       DATA ( BM3I(  8, 5,IBETA ), IBETA = 1, 10 ) /
     + 0.87576,0.88323,0.89935,0.91799,0.93583,0.95126,0.96377,0.97345,
     + 0.98072,0.98606/
       DATA ( BM3I(  8, 6,IBETA ), IBETA = 1, 10 ) /
     + 0.83078,0.83894,0.85705,0.87899,0.90126,0.92179,0.93950,0.95404,
     + 0.96551,0.97430/
       DATA ( BM3I(  8, 7,IBETA ), IBETA = 1, 10 ) /
     + 0.85727,0.86294,0.87558,0.89111,0.90723,0.92260,0.93645,0.94841,
     + 0.95838,0.96643/
       DATA ( BM3I(  8, 8,IBETA ), IBETA = 1, 10 ) /
     + 0.93337,0.93615,0.94220,0.94937,0.95647,0.96292,0.96840,0.97283,
     + 0.97619,0.97854/
       DATA ( BM3I(  8, 9,IBETA ), IBETA = 1, 10 ) /
     + 0.97790,0.97891,0.98105,0.98346,0.98569,0.98751,0.98879,0.98950,
     + 0.98961,0.98912/
       DATA ( BM3I(  8,10,IBETA ), IBETA = 1, 10 ) /
     + 0.99337,0.99367,0.99430,0.99493,0.99541,0.99562,0.99551,0.99501,
     + 0.99410,0.99274/
       DATA ( BM3I(  9, 1,IBETA ), IBETA = 1, 10 ) /
     + 0.98470,0.98594,0.98844,0.99106,0.99334,0.99514,0.99650,0.99749,
     + 0.99821,0.99872/
       DATA ( BM3I(  9, 2,IBETA ), IBETA = 1, 10 ) /
     + 0.98184,0.98330,0.98624,0.98934,0.99205,0.99420,0.99582,0.99701,
     + 0.99787,0.99848/
       DATA ( BM3I(  9, 3,IBETA ), IBETA = 1, 10 ) /
     + 0.97288,0.97498,0.97927,0.98385,0.98789,0.99113,0.99360,0.99541,
     + 0.99673,0.99766/
       DATA ( BM3I(  9, 4,IBETA ), IBETA = 1, 10 ) /
     + 0.95403,0.95741,0.96440,0.97202,0.97887,0.98444,0.98872,0.99190,
     + 0.99421,0.99586/
       DATA ( BM3I(  9, 5,IBETA ), IBETA = 1, 10 ) /
     + 0.91845,0.92399,0.93567,0.94873,0.96076,0.97079,0.97865,0.98457,
     + 0.98892,0.99206/
       DATA ( BM3I(  9, 6,IBETA ), IBETA = 1, 10 ) /
     + 0.86762,0.87533,0.89202,0.91148,0.93027,0.94669,0.96013,0.97062,
     + 0.97855,0.98441/
       DATA ( BM3I(  9, 7,IBETA ), IBETA = 1, 10 ) /
     + 0.84550,0.85253,0.86816,0.88721,0.90671,0.92490,0.94083,0.95413,
     + 0.96481,0.97314/
       DATA ( BM3I(  9, 8,IBETA ), IBETA = 1, 10 ) /
     + 0.90138,0.90544,0.91437,0.92513,0.93602,0.94615,0.95506,0.96258,
     + 0.96868,0.97347/
       DATA ( BM3I(  9, 9,IBETA ), IBETA = 1, 10 ) /
     + 0.96248,0.96415,0.96773,0.97187,0.97583,0.97925,0.98198,0.98394,
     + 0.98514,0.98559/
       DATA ( BM3I(  9,10,IBETA ), IBETA = 1, 10 ) /
     + 0.98837,0.98892,0.99005,0.99127,0.99232,0.99306,0.99339,0.99328,
     + 0.99269,0.99161/
       DATA ( BM3I( 10, 1,IBETA ), IBETA = 1, 10 ) /
     + 0.99080,0.99158,0.99311,0.99471,0.99607,0.99715,0.99795,0.99853,
     + 0.99895,0.99925/
       DATA ( BM3I( 10, 2,IBETA ), IBETA = 1, 10 ) /
     + 0.98910,0.99001,0.99182,0.99371,0.99533,0.99661,0.99757,0.99826,
     + 0.99876,0.99912/
       DATA ( BM3I( 10, 3,IBETA ), IBETA = 1, 10 ) /
     + 0.98374,0.98506,0.98772,0.99051,0.99294,0.99486,0.99630,0.99736,
     + 0.99812,0.99866/
       DATA ( BM3I( 10, 4,IBETA ), IBETA = 1, 10 ) /
     + 0.97238,0.97453,0.97892,0.98361,0.98773,0.99104,0.99354,0.99538,
     + 0.99671,0.99765/
       DATA ( BM3I( 10, 5,IBETA ), IBETA = 1, 10 ) /
     + 0.94961,0.95333,0.96103,0.96941,0.97693,0.98303,0.98772,0.99119,
     + 0.99371,0.99551/
       DATA ( BM3I( 10, 6,IBETA ), IBETA = 1, 10 ) /
     + 0.90943,0.91550,0.92834,0.94275,0.95608,0.96723,0.97600,0.98263,
     + 0.98751,0.99103/
       DATA ( BM3I( 10, 7,IBETA ), IBETA = 1, 10 ) /
     + 0.86454,0.87200,0.88829,0.90749,0.92630,0.94300,0.95687,0.96785,
     + 0.97626,0.98254/
       DATA ( BM3I( 10, 8,IBETA ), IBETA = 1, 10 ) /
     + 0.87498,0.88048,0.89264,0.90737,0.92240,0.93642,0.94877,0.95917,
     + 0.96762,0.97429/
       DATA ( BM3I( 10, 9,IBETA ), IBETA = 1, 10 ) /
     + 0.93946,0.94209,0.94781,0.95452,0.96111,0.96704,0.97203,0.97602,
     + 0.97900,0.98106/
       DATA ( BM3I( 10,10,IBETA ), IBETA = 1, 10 ) /
     + 0.97977,0.98071,0.98270,0.98492,0.98695,0.98858,0.98970,0.99027,
     + 0.99026,0.98968/

C FSB FM correction for INTRAmodal M2 coagulation  
       DATA BM2II /
     &  0.707107,  0.720583,  0.745310,  0.748056,  0.696935,
     &  0.604164,  0.504622,  0.416559,  0.343394,  0.283641/

C *** total correction for INTRAmodal M2 coagulation
      DATA BM2IITT /
     &  1.000000,  0.907452,  0.680931,  0.409815,  0.196425,
     &  0.078814,  0.028473,  0.009800,  0.003322,  0.001129/

C FSB FM correction for M2 i to j coagulation
      DATA ( BM2IJ(  1,  1,IBETA ), IBETA = 1, 10 ) /
     &  0.707107,  0.716828,  0.738240,  0.764827,  0.793610,
     &  0.822843,  0.851217,  0.877670,  0.901404,  0.921944/
      DATA ( BM2IJ(  1,  2,IBETA ), IBETA = 1, 10 ) /
     &  0.719180,  0.727975,  0.747638,  0.772334,  0.799234,
     &  0.826666,  0.853406,  0.878482,  0.901162,  0.920987/
      DATA ( BM2IJ(  1,  3,IBETA ), IBETA = 1, 10 ) /
     &  0.760947,  0.767874,  0.783692,  0.803890,  0.826015,
     &  0.848562,  0.870498,  0.891088,  0.909823,  0.926400/
      DATA ( BM2IJ(  1,  4,IBETA ), IBETA = 1, 10 ) /
     &  0.830926,  0.836034,  0.847708,  0.862528,  0.878521,
     &  0.894467,  0.909615,  0.923520,  0.935959,  0.946858/
      DATA ( BM2IJ(  1,  5,IBETA ), IBETA = 1, 10 ) /
     &  0.903643,  0.907035,  0.914641,  0.924017,  0.933795,
     &  0.943194,  0.951806,  0.959449,  0.966087,  0.971761/
      DATA ( BM2IJ(  1,  6,IBETA ), IBETA = 1, 10 ) /
     &  0.954216,  0.956094,  0.960211,  0.965123,  0.970068,
     &  0.974666,  0.978750,  0.982277,  0.985268,  0.987775/
      DATA ( BM2IJ(  1,  7,IBETA ), IBETA = 1, 10 ) /
     &  0.980546,  0.981433,  0.983343,  0.985568,  0.987751,
     &  0.989735,  0.991461,  0.992926,  0.994150,  0.995164/
      DATA ( BM2IJ(  1,  8,IBETA ), IBETA = 1, 10 ) /
     &  0.992142,  0.992524,  0.993338,  0.994272,  0.995174,
     &  0.995981,  0.996675,  0.997257,  0.997740,  0.998137/
      DATA ( BM2IJ(  1,  9,IBETA ), IBETA = 1, 10 ) /
     &  0.996868,  0.997026,  0.997361,  0.997742,  0.998106,
     &  0.998430,  0.998705,  0.998935,  0.999125,  0.999280/
      DATA ( BM2IJ(  1, 10,IBETA ), IBETA = 1, 10 ) /
     &  0.998737,  0.998802,  0.998939,  0.999094,  0.999241,
     &  0.999371,  0.999481,  0.999573,  0.999648,  0.999709/
      DATA ( BM2IJ(  2,  1,IBETA ), IBETA = 1, 10 ) /
     &  0.729600,  0.739948,  0.763059,  0.791817,  0.822510,
     &  0.852795,  0.881000,  0.905999,  0.927206,  0.944532/
      DATA ( BM2IJ(  2,  2,IBETA ), IBETA = 1, 10 ) /
     &  0.727025,  0.737116,  0.759615,  0.787657,  0.817740,
     &  0.847656,  0.875801,  0.901038,  0.922715,  0.940643/
      DATA ( BM2IJ(  2,  3,IBETA ), IBETA = 1, 10 ) /
     &  0.738035,  0.746779,  0.766484,  0.791340,  0.818324,
     &  0.845546,  0.871629,  0.895554,  0.916649,  0.934597/
      DATA ( BM2IJ(  2,  4,IBETA ), IBETA = 1, 10 ) /
     &  0.784185,  0.790883,  0.806132,  0.825501,  0.846545,
     &  0.867745,  0.888085,  0.906881,  0.923705,  0.938349/
      DATA ( BM2IJ(  2,  5,IBETA ), IBETA = 1, 10 ) /
     &  0.857879,  0.862591,  0.873238,  0.886539,  0.900645,
     &  0.914463,  0.927360,  0.939004,  0.949261,  0.958125/
      DATA ( BM2IJ(  2,  6,IBETA ), IBETA = 1, 10 ) /
     &  0.925441,  0.928304,  0.934645,  0.942324,  0.950181,
     &  0.957600,  0.964285,  0.970133,  0.975147,  0.979388/
      DATA ( BM2IJ(  2,  7,IBETA ), IBETA = 1, 10 ) /
     &  0.966728,  0.968176,  0.971323,  0.975027,  0.978705,
     &  0.982080,  0.985044,  0.987578,  0.989710,  0.991485/
      DATA ( BM2IJ(  2,  8,IBETA ), IBETA = 1, 10 ) /
     &  0.986335,  0.986980,  0.988362,  0.989958,  0.991511,
     &  0.992912,  0.994122,  0.995143,  0.995992,  0.996693/
      DATA ( BM2IJ(  2,  9,IBETA ), IBETA = 1, 10 ) /
     &  0.994547,  0.994817,  0.995391,  0.996046,  0.996677,
     &  0.997238,  0.997719,  0.998122,  0.998454,  0.998727/
      DATA ( BM2IJ(  2, 10,IBETA ), IBETA = 1, 10 ) /
     &  0.997817,  0.997928,  0.998163,  0.998429,  0.998683,
     &  0.998908,  0.999099,  0.999258,  0.999389,  0.999497/
      DATA ( BM2IJ(  3,  1,IBETA ), IBETA = 1, 10 ) /
     &  0.783612,  0.793055,  0.814468,  0.841073,  0.868769,
     &  0.894963,  0.918118,  0.937527,  0.953121,  0.965244/
      DATA ( BM2IJ(  3,  2,IBETA ), IBETA = 1, 10 ) /
     &  0.772083,  0.781870,  0.803911,  0.831238,  0.859802,
     &  0.887036,  0.911349,  0.931941,  0.948649,  0.961751/
      DATA ( BM2IJ(  3,  3,IBETA ), IBETA = 1, 10 ) /
     &  0.755766,  0.765509,  0.787380,  0.814630,  0.843526,
     &  0.871670,  0.897443,  0.919870,  0.938557,  0.953576/
      DATA ( BM2IJ(  3,  4,IBETA ), IBETA = 1, 10 ) /
     &  0.763816,  0.772145,  0.790997,  0.814784,  0.840434,
     &  0.865978,  0.890034,  0.911671,  0.930366,  0.945963/
      DATA ( BM2IJ(  3,  5,IBETA ), IBETA = 1, 10 ) /
     &  0.813597,  0.819809,  0.833889,  0.851618,  0.870640,
     &  0.889514,  0.907326,  0.923510,  0.937768,  0.950003/
      DATA ( BM2IJ(  3,  6,IBETA ), IBETA = 1, 10 ) /
     &  0.886317,  0.890437,  0.899643,  0.910955,  0.922730,
     &  0.934048,  0.944422,  0.953632,  0.961624,  0.968444/
      DATA ( BM2IJ(  3,  7,IBETA ), IBETA = 1, 10 ) /
     &  0.944565,  0.946855,  0.951872,  0.957854,  0.963873,
     &  0.969468,  0.974438,  0.978731,  0.982372,  0.985424/
      DATA ( BM2IJ(  3,  8,IBETA ), IBETA = 1, 10 ) /
     &  0.976358,  0.977435,  0.979759,  0.982467,  0.985125,
     &  0.987540,  0.989642,  0.991425,  0.992916,  0.994150/
      DATA ( BM2IJ(  3,  9,IBETA ), IBETA = 1, 10 ) /
     &  0.990471,  0.990932,  0.991917,  0.993048,  0.994142,
     &  0.995121,  0.995964,  0.996671,  0.997258,  0.997740/
      DATA ( BM2IJ(  3, 10,IBETA ), IBETA = 1, 10 ) /
     &  0.996199,  0.996389,  0.996794,  0.997254,  0.997694,
     &  0.998086,  0.998420,  0.998699,  0.998929,  0.999117/
      DATA ( BM2IJ(  4,  1,IBETA ), IBETA = 1, 10 ) /
     &  0.844355,  0.852251,  0.869914,  0.891330,  0.912823,
     &  0.932259,  0.948642,  0.961767,  0.971897,  0.979510/
      DATA ( BM2IJ(  4,  2,IBETA ), IBETA = 1, 10 ) /
     &  0.831550,  0.839954,  0.858754,  0.881583,  0.904592,
     &  0.925533,  0.943309,  0.957647,  0.968779,  0.977185/
      DATA ( BM2IJ(  4,  3,IBETA ), IBETA = 1, 10 ) /
     &  0.803981,  0.813288,  0.834060,  0.859400,  0.885285,
     &  0.909286,  0.930084,  0.947193,  0.960714,  0.971078/
      DATA ( BM2IJ(  4,  4,IBETA ), IBETA = 1, 10 ) /
     &  0.781787,  0.791080,  0.811931,  0.837749,  0.864768,
     &  0.890603,  0.913761,  0.933477,  0.949567,  0.962261/
      DATA ( BM2IJ(  4,  5,IBETA ), IBETA = 1, 10 ) /
     &  0.791591,  0.799355,  0.816916,  0.838961,  0.862492,
     &  0.885595,  0.907003,  0.925942,  0.942052,  0.955310/
      DATA ( BM2IJ(  4,  6,IBETA ), IBETA = 1, 10 ) /
     &  0.844933,  0.850499,  0.863022,  0.878593,  0.895038,
     &  0.911072,  0.925939,  0.939227,  0.950765,  0.960550/
      DATA ( BM2IJ(  4,  7,IBETA ), IBETA = 1, 10 ) /
     &  0.912591,  0.916022,  0.923607,  0.932777,  0.942151,
     &  0.951001,  0.958976,  0.965950,  0.971924,  0.976965/
      DATA ( BM2IJ(  4,  8,IBETA ), IBETA = 1, 10 ) /
     &  0.959859,  0.961617,  0.965433,  0.969924,  0.974382,
     &  0.978472,  0.982063,  0.985134,  0.987716,  0.989865/
      DATA ( BM2IJ(  4,  9,IBETA ), IBETA = 1, 10 ) /
     &  0.983377,  0.984162,  0.985844,  0.987788,  0.989681,
     &  0.991386,  0.992860,  0.994104,  0.995139,  0.995991/
      DATA ( BM2IJ(  4, 10,IBETA ), IBETA = 1, 10 ) /
     &  0.993343,  0.993672,  0.994370,  0.995169,  0.995937,
     &  0.996622,  0.997209,  0.997700,  0.998106,  0.998439/
      DATA ( BM2IJ(  5,  1,IBETA ), IBETA = 1, 10 ) /
     &  0.895806,  0.901918,  0.915233,  0.930783,  0.945768,
     &  0.958781,  0.969347,  0.977540,  0.983697,  0.988225/
      DATA ( BM2IJ(  5,  2,IBETA ), IBETA = 1, 10 ) /
     &  0.885634,  0.892221,  0.906629,  0.923540,  0.939918,
     &  0.954213,  0.965873,  0.974951,  0.981794,  0.986840/
      DATA ( BM2IJ(  5,  3,IBETA ), IBETA = 1, 10 ) /
     &  0.860120,  0.867858,  0.884865,  0.904996,  0.924724,
     &  0.942177,  0.956602,  0.967966,  0.976616,  0.983043/
      DATA ( BM2IJ(  5,  4,IBETA ), IBETA = 1, 10 ) /
     &  0.827462,  0.836317,  0.855885,  0.879377,  0.902897,
     &  0.924232,  0.942318,  0.956900,  0.968222,  0.976774/
      DATA ( BM2IJ(  5,  5,IBETA ), IBETA = 1, 10 ) /
     &  0.805527,  0.814279,  0.833853,  0.857892,  0.882726,
     &  0.906095,  0.926690,  0.943938,  0.957808,  0.968615/
      DATA ( BM2IJ(  5,  6,IBETA ), IBETA = 1, 10 ) /
     &  0.820143,  0.827223,  0.843166,  0.863002,  0.883905,
     &  0.904128,  0.922585,  0.938687,  0.952222,  0.963255/
      DATA ( BM2IJ(  5,  7,IBETA ), IBETA = 1, 10 ) /
     &  0.875399,  0.880208,  0.890929,  0.904065,  0.917699,
     &  0.930756,  0.942656,  0.953131,  0.962113,  0.969657/
      DATA ( BM2IJ(  5,  8,IBETA ), IBETA = 1, 10 ) /
     &  0.934782,  0.937520,  0.943515,  0.950656,  0.957840,
     &  0.964516,  0.970446,  0.975566,  0.979905,  0.983534/
      DATA ( BM2IJ(  5,  9,IBETA ), IBETA = 1, 10 ) /
     &  0.971369,  0.972679,  0.975505,  0.978797,  0.982029,
     &  0.984964,  0.987518,  0.989685,  0.991496,  0.992994/
      DATA ( BM2IJ(  5, 10,IBETA ), IBETA = 1, 10 ) /
     &  0.988329,  0.988893,  0.990099,  0.991485,  0.992825,
     &  0.994025,  0.995058,  0.995925,  0.996643,  0.997234/
      DATA ( BM2IJ(  6,  1,IBETA ), IBETA = 1, 10 ) /
     &  0.933384,  0.937784,  0.947130,  0.957655,  0.967430,
     &  0.975639,  0.982119,  0.987031,  0.990657,  0.993288/
      DATA ( BM2IJ(  6,  2,IBETA ), IBETA = 1, 10 ) /
     &  0.926445,  0.931227,  0.941426,  0.952975,  0.963754,
     &  0.972845,  0.980044,  0.985514,  0.989558,  0.992498/
      DATA ( BM2IJ(  6,  3,IBETA ), IBETA = 1, 10 ) /
     &  0.907835,  0.913621,  0.926064,  0.940308,  0.953745,
     &  0.965189,  0.974327,  0.981316,  0.986510,  0.990297/
      DATA ( BM2IJ(  6,  4,IBETA ), IBETA = 1, 10 ) /
     &  0.879088,  0.886306,  0.901945,  0.920079,  0.937460,
     &  0.952509,  0.964711,  0.974166,  0.981265,  0.986484/
      DATA ( BM2IJ(  6,  5,IBETA ), IBETA = 1, 10 ) /
     &  0.846500,  0.854862,  0.873189,  0.894891,  0.916264,
     &  0.935315,  0.951197,  0.963812,  0.973484,  0.980715/
      DATA ( BM2IJ(  6,  6,IBETA ), IBETA = 1, 10 ) /
     &  0.828137,  0.836250,  0.854310,  0.876287,  0.898710,
     &  0.919518,  0.937603,  0.952560,  0.964461,  0.973656/
      DATA ( BM2IJ(  6,  7,IBETA ), IBETA = 1, 10 ) /
     &  0.848595,  0.854886,  0.868957,  0.886262,  0.904241,
     &  0.921376,  0.936799,  0.950096,  0.961172,  0.970145/
      DATA ( BM2IJ(  6,  8,IBETA ), IBETA = 1, 10 ) /
     &  0.902919,  0.906922,  0.915760,  0.926427,  0.937312,
     &  0.947561,  0.956758,  0.964747,  0.971525,  0.977175/
      DATA ( BM2IJ(  6,  9,IBETA ), IBETA = 1, 10 ) /
     &  0.952320,  0.954434,  0.959021,  0.964418,  0.969774,
     &  0.974688,  0.979003,  0.982690,  0.985789,  0.988364/
      DATA ( BM2IJ(  6, 10,IBETA ), IBETA = 1, 10 ) /
     &  0.979689,  0.980650,  0.982712,  0.985093,  0.987413,
     &  0.989502,  0.991308,  0.992831,  0.994098,  0.995142/
      DATA ( BM2IJ(  7,  1,IBETA ), IBETA = 1, 10 ) /
     &  0.958611,  0.961598,  0.967817,  0.974620,  0.980752,
     &  0.985771,  0.989650,  0.992543,  0.994653,  0.996171/
      DATA ( BM2IJ(  7,  2,IBETA ), IBETA = 1, 10 ) /
     &  0.954225,  0.957488,  0.964305,  0.971795,  0.978576,
     &  0.984144,  0.988458,  0.991681,  0.994034,  0.995728/
      DATA ( BM2IJ(  7,  3,IBETA ), IBETA = 1, 10 ) /
     &  0.942147,  0.946158,  0.954599,  0.963967,  0.972529,
     &  0.979612,  0.985131,  0.989271,  0.992301,  0.994487/
      DATA ( BM2IJ(  7,  4,IBETA ), IBETA = 1, 10 ) /
     &  0.921821,  0.927048,  0.938140,  0.950598,  0.962118,
     &  0.971752,  0.979326,  0.985046,  0.989254,  0.992299/
      DATA ( BM2IJ(  7,  5,IBETA ), IBETA = 1, 10 ) /
     &  0.893419,  0.900158,  0.914598,  0.931070,  0.946584,
     &  0.959795,  0.970350,  0.978427,  0.984432,  0.988811/
      DATA ( BM2IJ(  7,  6,IBETA ), IBETA = 1, 10 ) /
     &  0.863302,  0.871111,  0.888103,  0.907990,  0.927305,
     &  0.944279,  0.958245,  0.969211,  0.977540,  0.983720/
      DATA ( BM2IJ(  7,  7,IBETA ), IBETA = 1, 10 ) /
     &  0.850182,  0.857560,  0.873890,  0.893568,  0.913408,
     &  0.931591,  0.947216,  0.960014,  0.970121,  0.977886/
      DATA ( BM2IJ(  7,  8,IBETA ), IBETA = 1, 10 ) /
     &  0.875837,  0.881265,  0.893310,  0.907936,  0.922910,
     &  0.936977,  0.949480,  0.960154,  0.968985,  0.976111/
      DATA ( BM2IJ(  7,  9,IBETA ), IBETA = 1, 10 ) /
     &  0.926228,  0.929445,  0.936486,  0.944868,  0.953293,
     &  0.961108,  0.968028,  0.973973,  0.978974,  0.983118/
      DATA ( BM2IJ(  7, 10,IBETA ), IBETA = 1, 10 ) /
     &  0.965533,  0.967125,  0.970558,  0.974557,  0.978484,
     &  0.982050,  0.985153,  0.987785,  0.989982,  0.991798/
      DATA ( BM2IJ(  8,  1,IBETA ), IBETA = 1, 10 ) /
     &  0.974731,  0.976674,  0.980660,  0.984926,  0.988689,
     &  0.991710,  0.994009,  0.995703,  0.996929,  0.997805/
      DATA ( BM2IJ(  8,  2,IBETA ), IBETA = 1, 10 ) /
     &  0.972062,  0.974192,  0.978571,  0.983273,  0.987432,
     &  0.990780,  0.993333,  0.995218,  0.996581,  0.997557/
      DATA ( BM2IJ(  8,  3,IBETA ), IBETA = 1, 10 ) /
     &  0.964662,  0.967300,  0.972755,  0.978659,  0.983921,
     &  0.988181,  0.991444,  0.993859,  0.995610,  0.996863/
      DATA ( BM2IJ(  8,  4,IBETA ), IBETA = 1, 10 ) /
     &  0.951782,  0.955284,  0.962581,  0.970559,  0.977737,
     &  0.983593,  0.988103,  0.991454,  0.993889,  0.995635/
      DATA ( BM2IJ(  8,  5,IBETA ), IBETA = 1, 10 ) /
     &  0.931947,  0.936723,  0.946751,  0.957843,  0.967942,
     &  0.976267,  0.982734,  0.987571,  0.991102,  0.993642/
      DATA ( BM2IJ(  8,  6,IBETA ), IBETA = 1, 10 ) /
     &  0.905410,  0.911665,  0.924950,  0.939908,  0.953798,
     &  0.965469,  0.974684,  0.981669,  0.986821,  0.990556/
      DATA ( BM2IJ(  8,  7,IBETA ), IBETA = 1, 10 ) /
     &  0.878941,  0.886132,  0.901679,  0.919688,  0.936970,
     &  0.951980,  0.964199,  0.973709,  0.980881,  0.986174/
      DATA ( BM2IJ(  8,  8,IBETA ), IBETA = 1, 10 ) /
     &  0.871653,  0.878218,  0.892652,  0.909871,  0.927034,
     &  0.942592,  0.955836,  0.966604,  0.975065,  0.981545/
      DATA ( BM2IJ(  8,  9,IBETA ), IBETA = 1, 10 ) /
     &  0.900693,  0.905239,  0.915242,  0.927232,  0.939335,
     &  0.950555,  0.960420,  0.968774,  0.975651,  0.981188/
      DATA ( BM2IJ(  8, 10,IBETA ), IBETA = 1, 10 ) /
     &  0.944922,  0.947435,  0.952894,  0.959317,  0.965689,
     &  0.971529,  0.976645,  0.981001,  0.984641,  0.987642/
      DATA ( BM2IJ(  9,  1,IBETA ), IBETA = 1, 10 ) /
     &  0.984736,  0.985963,  0.988453,  0.991078,  0.993357,
     &  0.995161,  0.996519,  0.997512,  0.998226,  0.998734/
      DATA ( BM2IJ(  9,  2,IBETA ), IBETA = 1, 10 ) /
     &  0.983141,  0.984488,  0.987227,  0.990119,  0.992636,
     &  0.994632,  0.996137,  0.997238,  0.998030,  0.998595/
      DATA ( BM2IJ(  9,  3,IBETA ), IBETA = 1, 10 ) /
     &  0.978726,  0.980401,  0.983819,  0.987450,  0.990626,
     &  0.993157,  0.995071,  0.996475,  0.997486,  0.998206/
      DATA ( BM2IJ(  9,  4,IBETA ), IBETA = 1, 10 ) /
     &  0.970986,  0.973224,  0.977818,  0.982737,  0.987072,
     &  0.990546,  0.993184,  0.995124,  0.996523,  0.997521/
      DATA ( BM2IJ(  9,  5,IBETA ), IBETA = 1, 10 ) /
     &  0.958579,  0.961700,  0.968149,  0.975116,  0.981307,
     &  0.986301,  0.990112,  0.992923,  0.994954,  0.996404/
      DATA ( BM2IJ(  9,  6,IBETA ), IBETA = 1, 10 ) /
     &  0.940111,  0.944479,  0.953572,  0.963506,  0.972436,
     &  0.979714,  0.985313,  0.989468,  0.992483,  0.994641/
      DATA ( BM2IJ(  9,  7,IBETA ), IBETA = 1, 10 ) /
     &  0.916127,  0.921878,  0.934003,  0.947506,  0.959899,
     &  0.970199,  0.978255,  0.984314,  0.988755,  0.991960/
      DATA ( BM2IJ(  9,  8,IBETA ), IBETA = 1, 10 ) /
     &  0.893848,  0.900364,  0.914368,  0.930438,  0.945700,
     &  0.958824,  0.969416,  0.977603,  0.983746,  0.988262/
      DATA ( BM2IJ(  9,  9,IBETA ), IBETA = 1, 10 ) /
     &  0.892161,  0.897863,  0.910315,  0.925021,  0.939523,
     &  0.952544,  0.963544,  0.972442,  0.979411,  0.984742/
      DATA ( BM2IJ(  9, 10,IBETA ), IBETA = 1, 10 ) /
     &  0.922260,  0.925966,  0.934047,  0.943616,  0.953152,
     &  0.961893,  0.969506,  0.975912,  0.981167,  0.985394/
      DATA ( BM2IJ( 10,  1,IBETA ), IBETA = 1, 10 ) /
     &  0.990838,  0.991598,  0.993128,  0.994723,  0.996092,
     &  0.997167,  0.997969,  0.998552,  0.998969,  0.999265/
      DATA ( BM2IJ( 10,  2,IBETA ), IBETA = 1, 10 ) /
     &  0.989892,  0.990727,  0.992411,  0.994167,  0.995678,
     &  0.996864,  0.997751,  0.998396,  0.998858,  0.999186/
      DATA ( BM2IJ( 10,  3,IBETA ), IBETA = 1, 10 ) /
     &  0.987287,  0.988327,  0.990428,  0.992629,  0.994529,
     &  0.996026,  0.997148,  0.997965,  0.998551,  0.998967/
      DATA ( BM2IJ( 10,  4,IBETA ), IBETA = 1, 10 ) /
     &  0.982740,  0.984130,  0.986952,  0.989926,  0.992508,
     &  0.994551,  0.996087,  0.997208,  0.998012,  0.998584/
      DATA ( BM2IJ( 10,  5,IBETA ), IBETA = 1, 10 ) /
     &  0.975380,  0.977330,  0.981307,  0.985529,  0.989216,
     &  0.992147,  0.994358,  0.995975,  0.997136,  0.997961/
      DATA ( BM2IJ( 10,  6,IBETA ), IBETA = 1, 10 ) /
     &  0.963911,  0.966714,  0.972465,  0.978614,  0.984022,
     &  0.988346,  0.991620,  0.994020,  0.995747,  0.996974/
      DATA ( BM2IJ( 10,  7,IBETA ), IBETA = 1, 10 ) /
     &  0.947187,  0.951161,  0.959375,  0.968258,  0.976160,
     &  0.982540,  0.987409,  0.991000,  0.993592,  0.995441/
      DATA ( BM2IJ( 10,  8,IBETA ), IBETA = 1, 10 ) /
     &  0.926045,  0.931270,  0.942218,  0.954297,  0.965273,
     &  0.974311,  0.981326,  0.986569,  0.990394,  0.993143/
      DATA ( BM2IJ( 10,  9,IBETA ), IBETA = 1, 10 ) /
     &  0.908092,  0.913891,  0.926288,  0.940393,  0.953667,
     &  0.964987,  0.974061,  0.981038,  0.986253,  0.990078/
      DATA ( BM2IJ( 10, 10,IBETA ), IBETA = 1, 10 ) /
     &  0.911143,  0.915972,  0.926455,  0.938721,  0.950701,
     &  0.961370,  0.970329,  0.977549,  0.983197,  0.987518/

C FSB Total correction factor for M2 coagulation j from i

      DATA ( BM2JI(  1, 1,IBETA ), IBETA = 1, 10 ) /
     &  0.753466,  0.756888,  0.761008,  0.759432,  0.748675,
     &  0.726951,  0.693964,  0.650915,  0.600227,  0.545000/
      DATA ( BM2JI(  1, 2,IBETA ), IBETA = 1, 10 ) /
     &  0.824078,  0.828698,  0.835988,  0.838943,  0.833454,
     &  0.817148,  0.789149,  0.750088,  0.701887,  0.647308/
      DATA ( BM2JI(  1, 3,IBETA ), IBETA = 1, 10 ) /
     &  1.007389,  1.014362,  1.028151,  1.041011,  1.047939,
     &  1.045707,  1.032524,  1.007903,  0.972463,  0.927667/
      DATA ( BM2JI(  1, 4,IBETA ), IBETA = 1, 10 ) /
     &  1.246157,  1.255135,  1.274249,  1.295351,  1.313362,
     &  1.325187,  1.329136,  1.324491,  1.311164,  1.289459/
      DATA ( BM2JI(  1, 5,IBETA ), IBETA = 1, 10 ) /
     &  1.450823,  1.459551,  1.478182,  1.499143,  1.518224,
     &  1.533312,  1.543577,  1.548882,  1.549395,  1.545364/
      DATA ( BM2JI(  1, 6,IBETA ), IBETA = 1, 10 ) /
     &  1.575248,  1.581832,  1.595643,  1.610866,  1.624601,
     &  1.635690,  1.643913,  1.649470,  1.652688,  1.653878/
      DATA ( BM2JI(  1, 7,IBETA ), IBETA = 1, 10 ) /
     &  1.638426,  1.642626,  1.651293,  1.660641,  1.668926,
     &  1.675571,  1.680572,  1.684147,  1.686561,  1.688047/
      DATA ( BM2JI(  1, 8,IBETA ), IBETA = 1, 10 ) /
     &  1.669996,  1.672392,  1.677283,  1.682480,  1.687028,
     &  1.690651,  1.693384,  1.695372,  1.696776,  1.697734/
      DATA ( BM2JI(  1, 9,IBETA ), IBETA = 1, 10 ) /
     &  1.686148,  1.687419,  1.689993,  1.692704,  1.695057,
     &  1.696922,  1.698329,  1.699359,  1.700099,  1.700621/
      DATA ( BM2JI(  1,10,IBETA ), IBETA = 1, 10 ) /
     &  1.694364,  1.695010,  1.696313,  1.697676,  1.698853,
     &  1.699782,  1.700482,  1.700996,  1.701366,  1.701631/
      DATA ( BM2JI(  2, 1,IBETA ), IBETA = 1, 10 ) /
     &  0.783166,  0.779369,  0.768044,  0.747572,  0.716709,
     &  0.675422,  0.624981,  0.567811,  0.507057,  0.445975/
      DATA ( BM2JI(  2, 2,IBETA ), IBETA = 1, 10 ) /
     &  0.848390,  0.847100,  0.840874,  0.826065,  0.800296,
     &  0.762625,  0.713655,  0.655545,  0.591603,  0.525571/
      DATA ( BM2JI(  2, 3,IBETA ), IBETA = 1, 10 ) /
     &  1.039894,  1.043786,  1.049445,  1.049664,  1.039407,
     &  1.015322,  0.975983,  0.922180,  0.856713,  0.783634/
      DATA ( BM2JI(  2, 4,IBETA ), IBETA = 1, 10 ) /
     &  1.345995,  1.356064,  1.376947,  1.398304,  1.412685,
     &  1.414611,  1.400652,  1.369595,  1.322261,  1.260993/
      DATA ( BM2JI(  2, 5,IBETA ), IBETA = 1, 10 ) /
     &  1.675575,  1.689859,  1.720957,  1.756659,  1.788976,
     &  1.812679,  1.824773,  1.824024,  1.810412,  1.784630/
      DATA ( BM2JI(  2, 6,IBETA ), IBETA = 1, 10 ) /
     &  1.919835,  1.933483,  1.962973,  1.996810,  2.028377,
     &  2.054172,  2.072763,  2.083963,  2.088190,  2.086052/
      DATA ( BM2JI(  2, 7,IBETA ), IBETA = 1, 10 ) /
     &  2.064139,  2.074105,  2.095233,  2.118909,  2.140688,
     &  2.158661,  2.172373,  2.182087,  2.188330,  2.191650/
      DATA ( BM2JI(  2, 8,IBETA ), IBETA = 1, 10 ) /
     &  2.144871,  2.150990,  2.163748,  2.177731,  2.190364,
     &  2.200712,  2.208687,  2.214563,  2.218716,  2.221502/
      DATA ( BM2JI(  2, 9,IBETA ), IBETA = 1, 10 ) /
     &  2.189223,  2.192595,  2.199540,  2.207033,  2.213706,
     &  2.219125,  2.223297,  2.226403,  2.228660,  2.230265/
      DATA ( BM2JI(  2,10,IBETA ), IBETA = 1, 10 ) /
     &  2.212595,  2.214342,  2.217912,  2.221723,  2.225082,
     &  2.227791,  2.229869,  2.231417,  2.232551,  2.233372/
      DATA ( BM2JI(  3, 1,IBETA ), IBETA = 1, 10 ) /
     &  0.837870,  0.824476,  0.793119,  0.750739,  0.700950,
     &  0.646691,  0.590508,  0.534354,  0.479532,  0.426856/
      DATA ( BM2JI(  3, 2,IBETA ), IBETA = 1, 10 ) /
     &  0.896771,  0.885847,  0.859327,  0.821694,  0.775312,
     &  0.722402,  0.665196,  0.605731,  0.545742,  0.486687/
      DATA ( BM2JI(  3, 3,IBETA ), IBETA = 1, 10 ) /
     &  1.076089,  1.071727,  1.058845,  1.036171,  1.002539,
     &  0.957521,  0.901640,  0.836481,  0.764597,  0.689151/
      DATA ( BM2JI(  3, 4,IBETA ), IBETA = 1, 10 ) /
     &  1.409571,  1.415168,  1.425346,  1.432021,  1.428632,
     &  1.409696,  1.371485,  1.312958,  1.236092,  1.145293/
      DATA ( BM2JI(  3, 5,IBETA ), IBETA = 1, 10 ) /
     &  1.862757,  1.880031,  1.918394,  1.963456,  2.004070,
     &  2.030730,  2.036144,  2.016159,  1.970059,  1.900079/
      DATA ( BM2JI(  3, 6,IBETA ), IBETA = 1, 10 ) /
     &  2.289741,  2.313465,  2.366789,  2.431612,  2.495597,
     &  2.549838,  2.588523,  2.608665,  2.609488,  2.591662/
      DATA ( BM2JI(  3, 7,IBETA ), IBETA = 1, 10 ) /
     &  2.597157,  2.618731,  2.666255,  2.722597,  2.777531,
     &  2.825187,  2.862794,  2.889648,  2.906199,  2.913380/
      DATA ( BM2JI(  3, 8,IBETA ), IBETA = 1, 10 ) /
     &  2.797975,  2.813116,  2.845666,  2.882976,  2.918289,
     &  2.948461,  2.972524,  2.990687,  3.003664,  3.012284/
      DATA ( BM2JI(  3, 9,IBETA ), IBETA = 1, 10 ) /
     &  2.920832,  2.929843,  2.948848,  2.970057,  2.989632,
     &  3.006057,  3.019067,  3.028979,  3.036307,  3.041574/
      DATA ( BM2JI(  3,10,IBETA ), IBETA = 1, 10 ) /
     &  2.989627,  2.994491,  3.004620,  3.015720,  3.025789,
     &  3.034121,  3.040664,  3.045641,  3.049347,  3.052066/
      DATA ( BM2JI(  4, 1,IBETA ), IBETA = 1, 10 ) /
     &  0.893179,  0.870897,  0.820996,  0.759486,  0.695488,
     &  0.634582,  0.579818,  0.532143,  0.490927,  0.454618/
      DATA ( BM2JI(  4, 2,IBETA ), IBETA = 1, 10 ) /
     &  0.948355,  0.927427,  0.880215,  0.821146,  0.758524,
     &  0.697680,  0.641689,  0.591605,  0.546919,  0.506208/
      DATA ( BM2JI(  4, 3,IBETA ), IBETA = 1, 10 ) /
     &  1.109562,  1.093648,  1.056438,  1.007310,  0.951960,
     &  0.894453,  0.837364,  0.781742,  0.727415,  0.673614/
      DATA ( BM2JI(  4, 4,IBETA ), IBETA = 1, 10 ) /
     &  1.423321,  1.417557,  1.402442,  1.379079,  1.347687,
     &  1.308075,  1.259703,  1.201983,  1.134778,  1.058878/
      DATA ( BM2JI(  4, 5,IBETA ), IBETA = 1, 10 ) /
     &  1.933434,  1.944347,  1.968765,  1.997653,  2.023054,
     &  2.036554,  2.029949,  1.996982,  1.934982,  1.845473/
      DATA ( BM2JI(  4, 6,IBETA ), IBETA = 1, 10 ) /
     &  2.547772,  2.577105,  2.645918,  2.735407,  2.830691,
     &  2.917268,  2.981724,  3.013684,  3.007302,  2.961560/
      DATA ( BM2JI(  4, 7,IBETA ), IBETA = 1, 10 ) /
     &  3.101817,  3.139271,  3.225851,  3.336402,  3.453409,
     &  3.563116,  3.655406,  3.724014,  3.766113,  3.781394/
      DATA ( BM2JI(  4, 8,IBETA ), IBETA = 1, 10 ) /
     &  3.540920,  3.573780,  3.647439,  3.737365,  3.828468,
     &  3.911436,  3.981317,  4.036345,  4.076749,  4.103751/
      DATA ( BM2JI(  4, 9,IBETA ), IBETA = 1, 10 ) /
     &  3.856771,  3.879363,  3.928579,  3.986207,  4.042173,
     &  4.091411,  4.132041,  4.164052,  4.188343,  4.206118/
      DATA ( BM2JI(  4,10,IBETA ), IBETA = 1, 10 ) /
     &  4.053923,  4.067191,  4.095509,  4.127698,  4.158037,
     &  4.184055,  4.205135,  4.221592,  4.234115,  4.243463/
      DATA ( BM2JI(  5, 1,IBETA ), IBETA = 1, 10 ) /
     &  0.935846,  0.906814,  0.843358,  0.768710,  0.695885,
     &  0.631742,  0.579166,  0.538471,  0.508410,  0.486863/
      DATA ( BM2JI(  5, 2,IBETA ), IBETA = 1, 10 ) /
     &  0.988308,  0.959524,  0.896482,  0.821986,  0.748887,
     &  0.684168,  0.630908,  0.589516,  0.558676,  0.536056/
      DATA ( BM2JI(  5, 3,IBETA ), IBETA = 1, 10 ) /
     &  1.133795,  1.107139,  1.048168,  0.977258,  0.906341,
     &  0.842477,  0.789093,  0.746731,  0.713822,  0.687495/
      DATA ( BM2JI(  5, 4,IBETA ), IBETA = 1, 10 ) /
     &  1.405692,  1.385781,  1.340706,  1.284776,  1.227085,
     &  1.173532,  1.127008,  1.087509,  1.052712,  1.018960/
      DATA ( BM2JI(  5, 5,IBETA ), IBETA = 1, 10 ) /
     &  1.884992,  1.879859,  1.868463,  1.854995,  1.841946,
     &  1.829867,  1.816972,  1.799319,  1.771754,  1.729406/
      DATA ( BM2JI(  5, 6,IBETA ), IBETA = 1, 10 ) /
     &  2.592275,  2.612268,  2.661698,  2.731803,  2.815139,
     &  2.901659,  2.978389,  3.031259,  3.048045,  3.021122/
      DATA ( BM2JI(  5, 7,IBETA ), IBETA = 1, 10 ) /
     &  3.390321,  3.435519,  3.545615,  3.698419,  3.876958,
     &  4.062790,  4.236125,  4.378488,  4.475619,  4.519170/
      DATA ( BM2JI(  5, 8,IBETA ), IBETA = 1, 10 ) /
     &  4.161376,  4.216558,  4.346896,  4.519451,  4.711107,
     &  4.902416,  5.077701,  5.226048,  5.341423,  5.421764/
      DATA ( BM2JI(  5, 9,IBETA ), IBETA = 1, 10 ) /
     &  4.843961,  4.892035,  5.001492,  5.138515,  5.281684,
     &  5.416805,  5.535493,  5.634050,  5.712063,  5.770996/
      DATA ( BM2JI(  5,10,IBETA ), IBETA = 1, 10 ) /
     &  5.352093,  5.385119,  5.458056,  5.545311,  5.632162,
     &  5.710566,  5.777005,  5.830863,  5.873123,  5.905442/
      DATA ( BM2JI(  6, 1,IBETA ), IBETA = 1, 10 ) /
     &  0.964038,  0.930794,  0.859433,  0.777776,  0.700566,
     &  0.634671,  0.582396,  0.543656,  0.517284,  0.501694/
      DATA ( BM2JI(  6, 2,IBETA ), IBETA = 1, 10 ) /
     &  1.013416,  0.979685,  0.907197,  0.824135,  0.745552,
     &  0.678616,  0.625870,  0.587348,  0.561864,  0.547674/
      DATA ( BM2JI(  6, 3,IBETA ), IBETA = 1, 10 ) /
     &  1.145452,  1.111457,  1.038152,  0.953750,  0.873724,
     &  0.805955,  0.753621,  0.717052,  0.694920,  0.684910/
      DATA ( BM2JI(  6, 4,IBETA ), IBETA = 1, 10 ) /
     &  1.376547,  1.345004,  1.276415,  1.196704,  1.121091,
     &  1.058249,  1.012197,  0.983522,  0.970323,  0.968933/
      DATA ( BM2JI(  6, 5,IBETA ), IBETA = 1, 10 ) /
     &  1.778801,  1.755897,  1.706074,  1.649008,  1.597602,
     &  1.560087,  1.540365,  1.538205,  1.549738,  1.568333/
      DATA ( BM2JI(  6, 6,IBETA ), IBETA = 1, 10 ) /
     &  2.447603,  2.445172,  2.443762,  2.451842,  2.475877,
     &  2.519039,  2.580118,  2.653004,  2.727234,  2.789738/
      DATA ( BM2JI(  6, 7,IBETA ), IBETA = 1, 10 ) /
     &  3.368490,  3.399821,  3.481357,  3.606716,  3.772101,
     &  3.969416,  4.184167,  4.396163,  4.582502,  4.721838/
      DATA ( BM2JI(  6, 8,IBETA ), IBETA = 1, 10 ) /
     &  4.426458,  4.489861,  4.648250,  4.877510,  5.160698,
     &  5.477495,  5.803123,  6.111250,  6.378153,  6.586050/
      DATA ( BM2JI(  6, 9,IBETA ), IBETA = 1, 10 ) /
     &  5.568061,  5.644988,  5.829837,  6.081532,  6.371214,
     &  6.672902,  6.963737,  7.226172,  7.449199,  7.627886/
      DATA ( BM2JI(  6,10,IBETA ), IBETA = 1, 10 ) /
     &  6.639152,  6.707020,  6.863974,  7.065285,  7.281744,
     &  7.492437,  7.683587,  7.847917,  7.983296,  8.090977/
      DATA ( BM2JI(  7, 1,IBETA ), IBETA = 1, 10 ) /
     &  0.980853,  0.945724,  0.871244,  0.787311,  0.708818,
     &  0.641987,  0.588462,  0.547823,  0.518976,  0.500801/
      DATA ( BM2JI(  7, 2,IBETA ), IBETA = 1, 10 ) /
     &  1.026738,  0.990726,  0.914306,  0.828140,  0.747637,
     &  0.679351,  0.625127,  0.584662,  0.556910,  0.540749/
      DATA ( BM2JI(  7, 3,IBETA ), IBETA = 1, 10 ) /
     &  1.146496,  1.108808,  1.028695,  0.938291,  0.854101,
     &  0.783521,  0.728985,  0.690539,  0.667272,  0.657977/
      DATA ( BM2JI(  7, 4,IBETA ), IBETA = 1, 10 ) /
     &  1.344846,  1.306434,  1.224543,  1.132031,  1.046571,
     &  0.976882,  0.926488,  0.896067,  0.884808,  0.891027/
      DATA ( BM2JI(  7, 5,IBETA ), IBETA = 1, 10 ) /
     &  1.670227,  1.634583,  1.558421,  1.472939,  1.396496,
     &  1.339523,  1.307151,  1.300882,  1.319622,  1.360166/
      DATA ( BM2JI(  7, 6,IBETA ), IBETA = 1, 10 ) /
     &  2.224548,  2.199698,  2.148284,  2.095736,  2.059319,
     &  2.050496,  2.075654,  2.136382,  2.229641,  2.347958/
      DATA ( BM2JI(  7, 7,IBETA ), IBETA = 1, 10 ) /
     &  3.104483,  3.105947,  3.118398,  3.155809,  3.230427,
     &  3.350585,  3.519071,  3.731744,  3.976847,  4.235616/
      DATA ( BM2JI(  7, 8,IBETA ), IBETA = 1, 10 ) /
     &  4.288426,  4.331456,  4.447024,  4.633023,  4.891991,
     &  5.221458,  5.610060,  6.036467,  6.471113,  6.880462/
      DATA ( BM2JI(  7, 9,IBETA ), IBETA = 1, 10 ) /
     &  5.753934,  5.837061,  6.048530,  6.363800,  6.768061,
     &  7.241280,  7.755346,  8.276666,  8.771411,  9.210826/
      DATA ( BM2JI(  7,10,IBETA ), IBETA = 1, 10 ) /
     &  7.466219,  7.568810,  7.819032,  8.168340,  8.582973,
     &  9.030174,  9.478159,  9.899834, 10.275940, 10.595910/
      DATA ( BM2JI(  8, 1,IBETA ), IBETA = 1, 10 ) /
     &  0.990036,  0.954782,  0.880531,  0.797334,  0.719410,
     &  0.652220,  0.596923,  0.552910,  0.519101,  0.494529/
      DATA ( BM2JI(  8, 2,IBETA ), IBETA = 1, 10 ) /
     &  1.032428,  0.996125,  0.919613,  0.833853,  0.753611,
     &  0.684644,  0.628260,  0.583924,  0.550611,  0.527407/
      DATA ( BM2JI(  8, 3,IBETA ), IBETA = 1, 10 ) /
     &  1.141145,  1.102521,  1.021017,  0.929667,  0.844515,
     &  0.772075,  0.714086,  0.670280,  0.639824,  0.621970/
      DATA ( BM2JI(  8, 4,IBETA ), IBETA = 1, 10 ) /
     &  1.314164,  1.273087,  1.186318,  1.089208,  0.999476,
     &  0.924856,  0.867948,  0.829085,  0.807854,  0.803759/
      DATA ( BM2JI(  8, 5,IBETA ), IBETA = 1, 10 ) /
     &  1.580611,  1.538518,  1.449529,  1.350459,  1.260910,
     &  1.190526,  1.143502,  1.121328,  1.124274,  1.151974/
      DATA ( BM2JI(  8, 6,IBETA ), IBETA = 1, 10 ) /
     &  2.016773,  1.977721,  1.895727,  1.806974,  1.732891,
     &  1.685937,  1.673026,  1.697656,  1.761039,  1.862391/
      DATA ( BM2JI(  8, 7,IBETA ), IBETA = 1, 10 ) /
     &  2.750093,  2.723940,  2.672854,  2.628264,  2.612250,
     &  2.640406,  2.723211,  2.866599,  3.071893,  3.335217/
      DATA ( BM2JI(  8, 8,IBETA ), IBETA = 1, 10 ) /
     &  3.881905,  3.887143,  3.913667,  3.981912,  4.111099,
     &  4.316575,  4.608146,  4.988157,  5.449592,  5.974848/
      DATA ( BM2JI(  8, 9,IBETA ), IBETA = 1, 10 ) /
     &  5.438870,  5.492742,  5.640910,  5.886999,  6.241641,
     &  6.710609,  7.289480,  7.960725,  8.693495,  9.446644/
      DATA ( BM2JI(  8,10,IBETA ), IBETA = 1, 10 ) /
     &  7.521152,  7.624621,  7.892039,  8.300444,  8.839787,
     &  9.493227, 10.231770, 11.015642, 11.799990, 12.542260/
      DATA ( BM2JI(  9, 1,IBETA ), IBETA = 1, 10 ) /
     &  0.994285,  0.960012,  0.887939,  0.807040,  0.730578,
     &  0.663410,  0.606466,  0.559137,  0.520426,  0.489429/
      DATA ( BM2JI(  9, 2,IBETA ), IBETA = 1, 10 ) /
     &  1.033505,  0.998153,  0.923772,  0.840261,  0.761383,
     &  0.692242,  0.633873,  0.585709,  0.546777,  0.516215/
      DATA ( BM2JI(  9, 3,IBETA ), IBETA = 1, 10 ) /
     &  1.132774,  1.094907,  1.015161,  0.925627,  0.841293,
     &  0.767888,  0.706741,  0.657439,  0.619135,  0.591119/
      DATA ( BM2JI(  9, 4,IBETA ), IBETA = 1, 10 ) /
     &  1.286308,  1.245273,  1.158809,  1.061889,  0.971208,
     &  0.893476,  0.830599,  0.782561,  0.748870,  0.729198/
      DATA ( BM2JI(  9, 5,IBETA ), IBETA = 1, 10 ) /
     &  1.511105,  1.467141,  1.374520,  1.271162,  1.175871,
     &  1.096887,  1.037243,  0.997820,  0.978924,  0.980962/
      DATA ( BM2JI(  9, 6,IBETA ), IBETA = 1, 10 ) /
     &  1.857468,  1.812177,  1.717002,  1.612197,  1.519171,
     &  1.448660,  1.405871,  1.393541,  1.413549,  1.467532/
      DATA ( BM2JI(  9, 7,IBETA ), IBETA = 1, 10 ) /
     &  2.430619,  2.388452,  2.301326,  2.210241,  2.139724,
     &  2.104571,  2.114085,  2.174696,  2.291294,  2.467500/
      DATA ( BM2JI(  9, 8,IBETA ), IBETA = 1, 10 ) /
     &  3.385332,  3.357690,  3.306611,  3.269804,  3.274462,
     &  3.340862,  3.484609,  3.717740,  4.048748,  4.481588/
      DATA ( BM2JI(  9, 9,IBETA ), IBETA = 1, 10 ) /
     &  4.850497,  4.858280,  4.896008,  4.991467,  5.171511,
     &  5.459421,  5.873700,  6.426128,  7.119061,  7.942603/
      DATA ( BM2JI(  9,10,IBETA ), IBETA = 1, 10 ) /
     &  6.957098,  7.020164,  7.197272,  7.499331,  7.946554,
     &  8.555048,  9.330503, 10.263610, 11.327454, 12.478332/
      DATA ( BM2JI( 10, 1,IBETA ), IBETA = 1, 10 ) /
     &  0.994567,  0.961842,  0.892854,  0.814874,  0.740198,
     &  0.673303,  0.615105,  0.565139,  0.522558,  0.486556/
      DATA ( BM2JI( 10, 2,IBETA ), IBETA = 1, 10 ) /
     &  1.031058,  0.997292,  0.926082,  0.845571,  0.768501,
     &  0.699549,  0.639710,  0.588538,  0.545197,  0.508894/
      DATA ( BM2JI( 10, 3,IBETA ), IBETA = 1, 10 ) /
     &  1.122535,  1.086287,  1.009790,  0.923292,  0.840626,
     &  0.766982,  0.703562,  0.650004,  0.605525,  0.569411/
      DATA ( BM2JI( 10, 4,IBETA ), IBETA = 1, 10 ) /
     &  1.261142,  1.221555,  1.137979,  1.043576,  0.953745,
     &  0.874456,  0.807292,  0.752109,  0.708326,  0.675477/
      DATA ( BM2JI( 10, 5,IBETA ), IBETA = 1, 10 ) /
     &  1.456711,  1.413432,  1.322096,  1.219264,  1.122319,
     &  1.038381,  0.969743,  0.916811,  0.879544,  0.858099/
      DATA ( BM2JI( 10, 6,IBETA ), IBETA = 1, 10 ) /
     &  1.741792,  1.695157,  1.596897,  1.487124,  1.385734,
     &  1.301670,  1.238638,  1.198284,  1.181809,  1.190689/
      DATA ( BM2JI( 10, 7,IBETA ), IBETA = 1, 10 ) /
     &  2.190197,  2.141721,  2.040226,  1.929245,  1.832051,
     &  1.760702,  1.721723,  1.719436,  1.757705,  1.840677/
      DATA ( BM2JI( 10, 8,IBETA ), IBETA = 1, 10 ) /
     &  2.940764,  2.895085,  2.801873,  2.707112,  2.638603,
     &  2.613764,  2.644686,  2.741255,  2.912790,  3.168519/
      DATA ( BM2JI( 10, 9,IBETA ), IBETA = 1, 10 ) /
     &  4.186191,  4.155844,  4.101953,  4.069102,  4.089886,
     &  4.189530,  4.389145,  4.707528,  5.161567,  5.765283/
      DATA ( BM2JI( 10,10,IBETA ), IBETA = 1, 10 ) /
     &  6.119526,  6.127611,  6.171174,  6.286528,  6.508738,
     &  6.869521,  7.396912,  8.113749,  9.034683, 10.162190/
        
C *** end of data statements.       

C----------------------------------------------------------------------

      CONSTII = ABS( HALF * ( TWO ) ** TWO3RDS - ONE )
      SQRTTWO = SQRT( TWO )      
      DLGSQT2 = ONE / LOG( SQRTTWO )

      ESAT01   = EXP( 0.125D0 * XXLSGAT * XXLSGAT )
      ESAC01   = EXP( 0.125D0 * XXLSGAC * XXLSGAC )

      ESAT04  = ESAT01 ** 4
      ESAC04  = ESAC01 ** 4

      ESAT05  = ESAT04 * ESAT01
      ESAC05  = ESAC04 * ESAC01

      ESAT08  = ESAT04 * ESAT04
      ESAC08  = ESAC04 * ESAC04

      ESAT09  = ESAT08 * ESAT01
      ESAC09  = ESAC08 * ESAC01

      ESAT16  = ESAT08 * ESAT08
      ESAC16  = ESAC08 * ESAC08

      ESAT20  = ESAT16 * ESAT04
      ESAC20  = ESAC16 * ESAC04

      ESAT24  = ESAT20 * ESAT04
      ESAC24  = ESAC20 * ESAC04

      ESAT25  = ESAT20 * ESAT05
      ESAC25  = ESAC20 * ESAC05

      ESAT36  = ESAT20 * ESAT16
      ESAC36  = ESAC20 * ESAC16

      ESAT49  = ESAT24 * ESAT25

      ESAT64  = ESAT20 * ESAT20 * ESAT24
      ESAC64  = ESAC20 * ESAC20 * ESAC24

      ESAT100 = ESAT64 * ESAT36

      DGAT2   = DGATK * DGATK
      DGAT3   = DGATK * DGATK * DGATK
      DGAC2   = DGACC * DGACC
      DGAC3   = DGACC * DGACC * DGACC

      SQDGAT  = SQRT( DGATK )
      SQDGAC  = SQRT( DGACC )
      SQDGAT5 = DGAT2 * SQDGAT
      SQDGAC5 = DGAC2 * SQDGAC
      SQDGAT7 = DGAT3 * SQDGAT

      XM2AT = DGAT2 * ESAT16
      XM3AT = DGAT3 * ESAT36

      XM2AC = DGAC2 * ESAC16
      XM3AC = DGAC3 * ESAC36

C *** For the free molecular regime:  Page H.3 of Whitby et al. (1991)
      R       = SQDGAC / SQDGAT
      R2      = R * R
      R3      = R2 * R
      R4      = R2 * R2
      R5      = R3 * R2
      R6      = R3 * R3
      R8      = R4 * R4
      RI1     = ONE / R
      RI2     = ONE / R2
      RI3     = ONE / R3
      RI4     = RI2 * RI2
      KNGAT   = TWO * LAMDA / DGATK
      KNGAC   = TWO * LAMDA / DGACC

C *** Calculate ratio of geometric mean diameters
      RAT = DGACC / DGATK

C *** Trap subscripts for BM0 and BM0I, between 1 and 10            
c     See page H.5 of Whitby et al. (1991)
      N2N = MAX( 1, MIN( 10, 
     &      NINT( 4.0 * ( SGATK - 0.75D0 ) ) ) )
      N2A = MAX( 1, MIN( 10, 
     &      NINT( 4.0 * ( SGACC - 0.75D0 ) ) ) )
      N1  = MAX( 1, MIN( 10,
     &       1 + NINT( DLGSQT2 * LOG( RAT ) ) ) )

C *** Intermodal Coagulation

C *** Set up for zeroth moment

C *** Near-continuum form:  Equation H.10a of Whitby et al. (1991)
         COAGNC0 = KNC * (         
     &    TWO + A * ( KNGAT * ( ESAT04 + R2 * ESAT16 * ESAC04 )
     &              + KNGAC * ( ESAC04 + RI2 * ESAC16 * ESAT04 ) )
     &              + ( R2 + RI2 ) * ESAT04 * ESAC04  ) 

     
C *** Free-molecular form:  Equation H.7a of Whitby et al. (1991)
         COAGFM0 = KFMATAC * SQDGAT * BM0IJ( N1,N2N,N2A ) * ( 
     &             ESAT01 + R * ESAC01 + TWO * R2 * ESAT01 * ESAC04 
     &           + R4 * ESAT09 * ESAC16 + RI3 * ESAT16 * ESAC09 
     &           + TWO * RI1 * ESAT04 + ESAC01  )  

C *** Loss to accumulation mode

C *** Harmonic Mean
      COAGATAC0 = COAGNC0 * COAGFM0 / ( COAGNC0 + COAGFM0 )
      
      QN12 = COAGATAC0

C *** Set up for second moment
C     The second moment equations are new and begin with equations A1
C     through A4 of Binkowski and Shankar (1995). After some algebraic
C     rearrangement and application of the extended mean value theorem
C     of integral calculus, equations are obtained that can be solved
C     analytically with correction factors as has been done by
C     Whitby et al. (1991)

C *** The term ( dp1 + dp2 ) ** (2/3) in Equations A3 and A4 of
C     Binkowski and Shankar (1995) is approximated by
C     (DGAT ** 3 + DGAC **3 ) ** 2/3 

C *** Near-continuum form
      I1NC = KNC * DGAT2 * (
     &       TWO * ESAT16
     &     + R2 * ESAT04 * ESAC04
     &     + RI2 * ESAT36 * ESAC04
     &     + A * KNGAT * (
     &           ESAT04
     &     +     RI2 * ESAT16 * ESAC04
     &     +     RI4 * ESAT36 * ESAC16
     &     +     R2 * ESAC04 )  )

C *** Free-molecular form
       I1FM = KFMATAC * SQDGAT5 * BM2IJ( N1,N2N,N2A ) * (
     &        ESAT25
     &      + TWO * R2 * ESAT09 * ESAC04
     &      + R4 * ESAT01 * ESAC16
     &      + RI3 * ESAT64 * ESAC09
     &      + TWO * RI1 * ESAT36 * ESAC01
     &      + R * ESAT16 * ESAC01  )

C *** Loss to accumulation mode

C *** Harmonic mean
      I1 = ( I1FM * I1NC ) / ( I1FM + I1NC )
      COAGATAC2 = I1 
      QS12 = COAGATAC2 

C *** Gain by accumulation mode
      COAGACAT2 = ( ( ONE + R6 ) ** TWO3RDS - R4 ) * I1 
      QS21 = COAGACAT2 * BM2JI( N1,N2N,N2A )

C *** Set up for third moment
     
C *** Near-continuum form: Equation H.10b of Whitby et al. (1991)
      COAGNC3 = KNC * DGAT3 * ( 
     &          TWO * ESAT36       
     &        + A * KNGAT * ( ESAT16 + R2 * ESAT04 * ESAC04 )                           
     &        + A * KNGAC * ( ESAT36 * ESAC04 + RI2 * ESAT64 * ESAC16 )
     &        + R2 * ESAT16 * ESAC04 + RI2 * ESAT64 * ESAC04 )
           
C *** Free_molecular form: Equation H.7b of Whitby et al. (1991)
      COAGFM3 = KFMATAC * SQDGAT7 * BM3I( N1, N2N, N2A ) * ( 
     &          ESAT49 
     &        + R * ESAT36  * ESAC01 
     &        + TWO * R2 * ESAT25  * ESAC04
     &        + R4 * ESAT09  * ESAC16
     &        + RI3 * ESAT100 * ESAC09
     &        + TWO * RI1 * ESAT64  * ESAC01 )   

C *** gain by accumulation mode = loss from Aitken mode

C *** Harmonic mean 
      COAGATAC3 = COAGNC3 * COAGFM3 / ( COAGNC3 + COAGFM3 )
      QV12 = COAGATAC3         

C *** Intramodal coagulation

C *** Zeroeth moment

C *** Aitken mode

C *** Near-continuum form: Equation H.12a of Whitby et al. (1991)
      COAGNC_AT = KNC * (ONE + ESAT08 + A * KNGAT * ( ESAT20 + ESAT04 ) )

C *** Free-molecular form: Equation H.11a of Whitby et al. (1991)
      COAGFM_AT = KFMAT * SQDGAT * BM0( N2N ) * 
     &           ( ESAT01 + ESAT25 + TWO * ESAT05 )

C *** Harmonic mean 
      COAGATAT0 = COAGFM_AT * COAGNC_AT / ( COAGFM_AT + COAGNC_AT )
      QN11 = COAGATAT0 
      
C *** Accumulation mode 

C *** Near-continuum form: Equation H.12a of Whitby et al. (1991)
      COAGNC_AC = KNC * (ONE + ESAC08 + A * KNGAC * ( ESAC20 + ESAC04 ) )

C *** Free-molecular form: Equation H.11a of Whitby et al. (1991)
      COAGFM_AC = KFMAC * SQDGAC * BM0( N2A ) * 
     &           ( ESAC01 + ESAC25 + TWO * ESAC05 ) 
       
C *** Harmonic mean 
      COAGACAC0 = COAGFM_AC * COAGNC_AC / ( COAGFM_AC + COAGNC_AC )
      QN22 = COAGACAC0

C *** Set up for second moment
C     The second moment equations are new and begin with 3.11a on Page
c     45 of Whitby et al. (1991). After some algebraic rearrangement and
c     application of the extended mean value theorem of integral calculus
c     equations are obtained that can be solved analytically with 
c     correction factors as has been done by Whitby et al. (1991)
           
C *** Aitken mode

C *** Near-continuum
      I1NC_AT = KNC * DGAT2 * (
     &          TWO * ESAT16
     &        + ESAT04 * ESAT04
     &        + ESAT36 * ESAT04
     &        + A * KNGAT * (
     &              TWO * ESAT04
     &        +     ESAT16 * ESAT04
     &        +     ESAT36 * ESAT16 )  )
     
C *** Free- molecular form
       I1FM_AT = KFMAT * SQDGAT5 * BM2II( N2N ) * (
     &           ESAT25
     &         +  TWO * ESAT09 * ESAT04
     &         +  ESAT01 * ESAT16
     &         +  ESAT64 * ESAT09
     &         +  TWO * ESAT36 * ESAT01
     &         +  ESAT16 * ESAT01  ) 

      I1_AT = ( I1NC_AT * I1FM_AT ) / ( I1NC_AT + I1FM_AT  )
      COAGATAT2 = CONSTII * I1_AT
      QS11 = COAGATAT2 * BM2IITT( N2N )

C *** Accumulation mode

C *** Near-continuum
      I1NC_AC = KNC * DGAC2 * (
     &          TWO * ESAC16
     &        + ESAC04 * ESAC04
     &        + ESAC36 * ESAC04
     &        + A * KNGAC * (
     &             TWO * ESAC04
     &        +    ESAC16 * ESAC04
     &        +    ESAC36 * ESAC16 )  )

C *** Free-molecular form
       I1FM_AC = KFMAC * SQDGAC5 * BM2II( N2A ) * (
     &           ESAC25
     &         + TWO * ESAC09 * ESAC04
     &         + ESAC01 * ESAC16
     &         + ESAC64 * ESAC09
     &         + TWO * ESAC36 * ESAC01
     &         + ESAC16 * ESAC01 )
          
      I1_AC = ( I1NC_AC * I1FM_AC ) / ( I1NC_AC + I1FM_AC  )
      COAGACAC2 = CONSTII * I1_AC
      QS22 = COAGACAC2 * BM2IITT( N2A )
     
      RETURN

      END SUBROUTINE GETCOAGS

C:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
      subroutine intracoag_gh( lamda, kfm, knc, dg, xlnsig,
     &                         quads11, quadn11 )

c  Gauss-Hermite numerical quadrature to calculate intramodal coagulation
c  rates for number and second moment
 
c FSB: This version runs in real*8 arithmetic
 
c *** This version calculates the coagulation coefficients using the harmonic
c     mean approach for both fm and nc cases.
c *** Does gauss-hermite quadrature for intra-modal coagulation integrals for
c     2nd moment for a lognormal distribution defined by dg, xlnsig,
c *** dg and xlnsig are the geometric mean diameters (meters) and the logarithms
c     of the geometric standard deviations (dimensionless) whose meaning is
c     defined below at the end of the routine ghxi, ghwi are the gauss-hermite
c     weights and n is one-half the number of abscissas, since an even number
c     of abscissas is used

      implicit none

c *** arguments:
      real( 8 ), intent( in )  :: lamda   ! mean free path
      real( 8 ), intent( in )  :: kfm, knc
      real,      intent( in )  :: dg, xlnsig
      real( 8 ), intent( out ) :: quads11, quadn11

c *** parameters:
      real( 8 ), parameter :: pi = 3.14159265358979324d0
      real( 8 ), parameter :: two3rds = 2.0d0 / 3.0d0
      real( 8 ), parameter :: sqrt2 = 1.414213562373095d0
      real( 8 ), parameter :: A = 1.246d0    ! approx Cunningham corr. factor
      real( 8 ), parameter :: twoA = 2.0d0 * A

c *** local variables:
      real( 8 ) :: sum1sfm, sum2sfm, sum1nfm, sum2nfm
      real( 8 ) :: sum1snc, sum2snc, sum1nnc, sum2nnc
      real( 8 ) :: xi, wxi, xf, dp1p, dp1m, dp1psq, dp1msq
      real( 8 ) :: v1p,v1m, a2p, a2m, v2p, v2m
      real( 8 ) :: yi, wyi, yf, dp2p, dp2m, dp2psq, dp2msq
      real( 8 ) :: dspp, dsmp, dspm, dsmm
      real( 8 ) :: bppfm, bmpfm, bpmfm, bmmfm
      real( 8 ) :: bppnc, bmpnc, bpmnc, bmmnc
      real( 8 ) :: xx1, xx2
      real( 8 ) :: xbsfm, xbsnc, xbnfm, xbnnc
      integer i, j

c *** Has a fixed number of Gauss-Herimite abscissas (n)
      integer, parameter :: n = 5   ! one-half the number of abscissas
      real( 8 ) :: ghxi( n ) ! Gauss-Hermite abscissas
      real( 8 ) :: ghwi( n ) ! Gauss-Hermite weights

c *** Values from Table 25.10 (page 924) of Abramowitz and Stegun,
c     Handbook of Mathematical Functions, National Bureau of Standards,
c     December 1965.
c    (Breaks in number to facilitate comparison with printed table.)

c *** tests show that 10 point is adquate.

      data ghxi/ 0.34290 13272 23705d0,
     &           1.03661 08297 89514d0,
     &           1.75668 36492 99882d0,
     &           2.53273 16742 32790d0,
     &           3.43615 91188 37738d0 /,

     &     ghwi/ 6.10862 63373 53d-01,
     &           2.40138 61108 23d-01,
     &           3.38743 94455 48d-02,
     &           1.34364 57467 81d-03,
     &           7.64043 28552 33d-06 /

c *** The following Statement Functions are based on expressions from
c     Binkowski & Shanker, Jour. Geophys. Research. Vol. 100, no. d12,
c     pp 26,191-26,209 December 20, 1995

      real( 8 ) :: betafm, betanc

c *** for Free Molecular, Eq. A5
      betafm( xx1, xx2 ) = kfm
     &                   * sqrt( 1.d0 / xx1 ** 3  + 1.d0 / xx2 ** 3 )
     &                   * ( xx1 + xx2 ) ** 2

c *** for Near Continuum, Eq. A6
      betanc( xx1, xx2 ) = knc * ( xx1 + xx2 )
     &                   * ( 1.0D0 / xx1 + 1.0d0 / xx2 + twoA * lamda
     &                      * ( 1.0d0 / xx1 ** 2 + 1.0d0 / xx2 ** 2 ) )
C-----------------------------------------------------------------------

      sum1sfm = 0.d0
      sum1snc = 0.d0
      sum1nfm = 0.d0
      sum1nnc = 0.d0
      do 201 i = 1, n

         sum2sfm = 0.d0
         sum2snc = 0.d0
         sum2nfm = 0.d0
         sum2nnc = 0.d0

         xi = ghxi( i )
         wxi = ghwi( i )
         xf = exp( sqrt2 * xi * xlnsig )
         dp1p = dg * xf
         dp1m = dg / xf
         dp1psq = dp1p * dp1p
         dp1msq = dp1m * dp1m
         v1p = dp1p * dp1psq
         v1m = dp1m * dp1msq
         do 101 j = 1, n
            yi = ghxi( j )
            wyi = ghwi( j )
            yf = exp( sqrt2 * yi * xlnsig )
            dp2p = dg * yf
            dp2m = dg / yf
            dp2psq = dp2p * dp2p
            dp2msq = dp2m * dp2m
            a2p = dp2psq
            a2m = dp2msq
            v2p = dp2p * dp2psq
            v2m = dp2m * dp2msq
            dspp = 0.5d0 * ( v1p + v2p ) ** two3rds - a2p
            dsmp = 0.5d0 * ( v1m + v2p ) ** two3rds - a2p
            dspm = 0.5d0 * ( v1p + v2m ) ** two3rds - a2m
            dsmm = 0.5d0 * ( v1m + v2m ) ** two3rds - a2m

            bppfm = betafm( dp1p, dp2p )
            bmpfm = betafm( dp1m, dp2p )
            bpmfm = betafm( dp1p, dp2m )
            bmmfm = betafm( dp1m, dp2m )

            bppnc = betanc( dp1p, dp2p )
            bmpnc = betanc( dp1m, dp2p )
            bpmnc = betanc( dp1p, dp2m )
            bmmnc = betanc( dp1m, dp2m )

            sum2sfm = sum2sfm + wyi * ( dspp * bppfm + dspm * bpmfm
     &                        + dsmp * bmpfm + dsmm * bmmfm )

            sum2nfm = sum2nfm + wyi * ( bppfm + bmpfm + bpmfm + bmmfm )
            sum2snc = sum2snc + wyi * ( dspp * bppnc + dspm * bpmnc
     &                        + dsmp * bmpnc + dsmm * bmmnc )
            sum2nnc = sum2nnc + wyi * ( bppnc + bmpnc + bpmnc + bmmnc )

101      continue
         sum1sfm = sum1sfm + wxi * sum2sfm
         sum1nfm = sum1nfm + wxi * sum2nfm

         sum1snc = sum1snc + wxi * sum2snc
         sum1nnc = sum1nnc + wxi * sum2nnc

201   continue

c *** quads11 is the intra-modal coagulation term for 2nd moment
      xbsfm = -sum1sfm / pi
      xbsnc = -sum1snc / pi
      quads11 = xbsfm * xbsnc / ( xbsfm + xbsnc )

c *** quadn11 is the intra-modal coagulation term for number
      xbnfm = 0.5D0 * sum1nfm / pi
      xbnnc = 0.5D0 * sum1nnc / pi
      quadn11 = xbnfm * xbnnc / ( xbnfm + xbnnc )

      return
      end subroutine intracoag_gh

C:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
      subroutine intercoag_gh( lamda, kfm, knc, dg1, dg2, xlnsig1, xlnsig2,
     &                         quads12, quads21, quadn12, quadv12 )

c  Gauss-Hermite numerical quadrature to calculate intermodal coagulation
c  for number, 2nd, and 3rd moment
 
c FSB This version runs in real*8 arithmetic
 
c *** This version calculates the coagulation coefficients using the harmonic
c     mean approach for both fm and nc cases.
c *** Does gauss-hermite quadrature for inter-modal coagulation integrals
c     for 2nd moment for two lognormal distributions defined by dg1, xlnsig1,
c     dg2, xlnsig2.
c *** dg and xlnsig are the geometric mean diameters (meters) and the logarithms
c     of the geometric standard deviations (dimensionless), whose meaning is
c     defined below at the end of the routine.
c     ghxi, ghwi are the gauss-hermite weights and n is one-half the
c     number of abscissas, since an even number of abscissas is used.

      implicit none


c *** arguments:
      real( 8 ), intent( in )  :: lamda   ! mean free path
      real( 8 ), intent( in )  :: kfm, knc
      real,      intent( in )  :: dg1, dg2, xlnsig1, xlnsig2
      real( 8 ), intent( out ) :: quads12, quads21, quadn12, quadv12

c *** parameters:
      real( 8 ), parameter :: pi = 3.14159265358979324d0
      real( 8 ), parameter :: two3rds = 2.0d0 / 3.0d0
      real( 8 ), parameter :: sqrt2 = 1.414213562373095d0
      real( 8 ), parameter :: A = 1.246d0    ! approx Cunningham corr. factor
      real( 8 ), parameter :: twoA = 2.0d0 * A

c *** local variables:
      real( 8 ) :: sum1s12fm, sum1s21fm, sum2s12fm, sum2s21fm
      real( 8 ) :: sum1nfm, sum2nfm
      real( 8 ) :: sum1vfm, sum2vfm
      real( 8 ) :: sum1s12nc, sum1s21nc, sum2s12nc, sum2s21nc
      real( 8 ) :: sum1nnc, sum2nnc
      real( 8 ) :: sum1vnc, sum2vnc
      real( 8 ) :: xi, wxi,xf, dp1p, dp1m, dp1psq, dp1msq
      real( 8 ) :: a1p, a1m, v1p, v1m
      real( 8 ) :: a2p, a2m, v2p, v2m
      real( 8 ) :: yi, wyi, yf, dp2p, dp2m, dp2psq, dp2msq
      real( 8 ) :: dspp, dsmp, dspm, dsmm
      real( 8 ) :: bppfm, bmpfm, bpmfm, bmmfm
      real( 8 ) :: bppnc, bmpnc, bpmnc, bmmnc
      real( 8 ) :: xx1, xx2
      real( 8 ) :: xbsfm, xbsnc, xbnfm, xbnnc, xbvfm, xbvnc
      integer i, j

c *** Has a fixed number of Gauss-Herimite abscissas (n)
      integer, parameter :: n = 5   ! one-half the number of abscissas
      real( 8 ) :: ghxi( n ) ! Gauss-Hermite abscissas
      real( 8 ) :: ghwi( n ) ! Gauss-Hermite weights

c *** Values from Table 25.10 (page 924) of Abramowitz and Stegun,
c     Handbook of Mathematical Functions, National Bureau of Standards,
c     December 1965.
c    (Breaks in number to facilitate comparison with printed table.)

c *** tests show that 10 point is adquate.

      data ghxi/ 0.34290 13272 23705D0,
     &           1.03661 08297 89514D0,
     &           1.75668 36492 99882D0,
     &           2.53273 16742 32790D0,
     &           3.43615 91188 37738D0 /,

     &     ghwi/ 6.10862 63373 53d-01,
     &           2.40138 61108 23d-01,
     &           3.38743 94455 48d-02,
     &           1.34364 57467 81d-03,
     &           7.64043 28552 33d-06 /

c *** The following Statement Functions are based on expressions from
c     Binkowski & Shanker, Jour. Geophys. Research. Vol. 100, no. d12,
c     pp 26,191-26,209 December 20, 1995

      real( 8 ) :: betafm, betanc

c *** for Free Molecular, Eq. A5
      betafm( xx1, xx2 ) = kfm
     &                   * sqrt( 1.d0 / xx1 ** 3  + 1.d0 / xx2 ** 3 )
     &                   * ( xx1 + xx2 ) ** 2

c *** for Near Continuum, Eq. A6
      betanc( xx1, xx2 ) = knc * ( xx1 + xx2 )
     &                   * ( 1.0D0 / xx1 + 1.0d0 / xx2 + twoA * lamda
     &                      * ( 1.0d0 / xx1 ** 2 + 1.0d0 / xx2 ** 2 ) )
c-----------------------------------------------------------------------

      sum1s12fm = 0.d0
      sum1s12nc = 0.d0
      sum1s21fm = 0.d0
      sum1s21nc = 0.d0
      sum1vnc = 0.d0
      sum1vfm = 0.d0
      sum1nfm = 0.d0
      sum1nnc = 0.d0
      do 201 i = 1, n

         sum2s12fm = 0.d0
         sum2s12nc = 0.d0
         sum2s21fm = 0.d0
         sum2s21nc = 0.d0
         sum2nfm = 0.d0
         sum2nnc = 0.d0
         sum2vnc = 0.d0
         sum2vfm = 0.d0
         xi = ghxi( i )
         wxi = ghwi( i )
         xf = exp( sqrt2 * xi *xlnsig1 )
         dp1p = dg1 * xf
         dp1m = dg1 / xf
         dp1psq = dp1p * dp1p
         dp1msq = dp1m * dp1m
         a1p = dp1psq
         a1m = dp1msq
         v1p = dp1p * dp1psq
         v1m = dp1m * dp1msq

         do 101 j = 1, n
            yi  = ghxi( j )
            wyi = ghwi( j )
            yf = exp( sqrt2 * yi * xlnsig2 )
            dp2p = dg2 * yf
            dp2m = dg2 / yf
            dp2psq = dp2p * dp2p
            dp2msq = dp2m * dp2m
            a2p  = dp2psq
            a2m  = dp2msq
            v2p  =  dp2p * dp2psq
            v2m  = dp2m * dp2msq
            dspp = ( v1p + v2p ) ** two3rds - a2p
            dsmp = ( v1m + v2p ) ** two3rds - a2p
            dspm = ( v1p + v2m ) ** two3rds - a2m
            dsmm = ( v1m + v2m ) ** two3rds - a2m

            bppfm = betafm( dp1p, dp2p )
            bmpfm = betafm( dp1m, dp2p )
            bpmfm = betafm( dp1p, dp2m )
            bmmfm = betafm( dp1m, dp2m )

            bppnc = betanc( dp1p, dp2p )
            bmpnc = betanc( dp1m, dp2p )
            bpmnc = betanc( dp1p, dp2m )
            bmmnc = betanc( dp1m, dp2m )

            sum2s12fm = sum2s12fm + wyi * ( a1p * bppfm + a1p * bpmfm
     &                            + a1m * bmpfm + a1m * bmmfm )

            sum2s21fm = sum2s21fm + wyi * ( dspp * bppfm + dspm * bpmfm
     &                            + dsmp * bmpfm + dsmm * bmmfm )

            sum2s12nc = sum2s12nc + wyi * ( a1p * bppnc + a1p * bpmnc
     &                            + a1m * bmpnc + a1m * bmmnc )
 
            sum2s21nc = sum2s21nc + wyi * ( dspp * bppnc + dspm * bpmnc
     &                            + dsmp * bmpnc + dsmm * bmmnc )

            sum2nfm = sum2nfm + wyi * ( bppfm + bmpfm + bpmfm + bmmfm )
  
            sum2nnc = sum2nnc + wyi * ( bppnc + bmpnc + bpmnc + bmmnc )

            sum2vfm = sum2vfm + wyi * ( v1p * ( bppfm + bpmfm )
     &                        + v1m * ( bmpfm + bmmfm ) )

            sum2vnc = sum2vnc + wyi * ( v1p * ( bppnc + bpmnc)
     &                        + v1m * ( bmpnc + bmmnc) )

101      continue

         sum1s12fm = sum1s12fm + wxi * sum2s12fm
         sum1s21fm = sum1s21fm + wxi * sum2s21fm
         sum1nfm   = sum1nfm + wxi * sum2nfm
         sum1vfm   = sum1vfm + wxi * sum2vfm

         sum1s12nc = sum1s12nc + wxi * sum2s12nc
         sum1s21nc = sum1s21nc + wxi * sum2s21nc
         sum1nnc   = sum1nnc + wxi * sum2nnc
         sum1vnc   = sum1vnc + wxi * sum2vnc

201   continue

C *** Second moment intermodal coagulation coefficients

c FSB NOTE: the transfer of second moment is not symmetric.
c     See equations A3 & A4 of Binkowski & Shankar (1995)

c *** to accumulation mode from Aitken mode
      xbsfm   = sum1s21fm / pi
      xbsnc   = sum1s21nc / pi
      quads21 = xbsfm * xbsnc / ( xbsfm + xbsnc )

c *** from Aitken mode to accumulation mode
      xbsfm   = sum1s12fm / pi
      xbsnc   = sum1s12nc / pi
      quads12 = xbsfm * xbsnc / ( xbsfm + xbsnc )

c *** quadn12 is the intermodal coagulation coefficient for number
      xbnfm   = sum1nfm / pi
      xbnnc   = sum1nnc / pi
      quadn12 = xbnfm * xbnnc / ( xbnfm + xbnnc )

c *** quadv12 is the intermodal coagulation coefficient for 3rd moment
      xbvfm = sum1vfm / pi
      xbvnc = sum1vnc / pi
      quadv12 = xbvfm * xbvnc / ( xbvfm + xbvnc )

      return
      end subroutine intercoag_gh

